#### **Tao Zhang**

Phase-field modeling of chemomechanical behavior in battery particles

PHASE-FIELD MODELING OF PHASE CHANGES AND MECHANICAL STRESSES IN ELECTRODE PARTICLES OF SECONDARY BATTERIES

SCHRIFTENREIHE DES INSTITUTS FÜR ANGEWANDTE MATERIALIEN

BAND 89

Tao Zhang

**Phase-field Modeling of Phase Changes and Mechanical Stresses in Electrode Particles of Secondary Batteries**

#### **Schriftenreihe des Instituts für Angewandte Materialien** *Band 89*

Karlsruher Institut für Technologie (KIT) Institut für Angewandte Materialien (IAM)

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buches.

# **Phase-field Modeling of Phase Changes and Mechanical Stresses in Electrode Particles of Secondary Batteries**

by Tao Zhang

Karlsruher Institut für Technologie Institut für Angewandte Materialien

Phase-field Modeling of Phase Changes and Mechanical Stresses in Electrode Particles of Secondary Batteries

Zur Erlangung des akademischen Grades eines Doktors der Ingenieurwissenschaften von der KIT-Fakultät für Maschinenbau des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Tao Zhang, M.Sc.

Tag der mündlichen Prüfung: 25. Juli 2019 Hauptreferent: Prof. Dr. Marc Kamlah Korreferent: Prof. Dr. Britta Nestler

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2021 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 2192-9963 ISBN 978-3-7315-1002-4 DOI 10.5445/KSP/1000100537

# **Abstract**

Electrochemical energy storage is needed for many mobile technical systems, such as communication and electromobility. Lithium-ion batteries (LIBs) have attracted intensive attention for electrochemical energy storage over the past decades. The enormous utilization of LIBs combined with the limited and unevenly distributed lithium source will drive up the price of lithium. In contrast to lithium, sodium-ion batteries (NIBs), which are based on the wide availability, abundance and low cost of sodium, become a potential promising alternative to LIBs. In most storage materials, phase changes occur during regular operation and, thus, cannot be avoided. The respective phases of a storage material possess different lattice constants giving rise to a strain mismatch of up to a few hundred percent which, in turn, causes mechanical stresses and, thus, leads to damage of the electrode particles. In view of this, it is the objective of this work to study phase changes and mechanical stresses in electrode particles by means of phase-field simulations. Thermodynamical phase-field models rely on continuous order parameters, thus, leading to diffuse interfaces between adjacent phases with no need for the cumbersome tracking of the position of a sharp interface.

In this work, for non-linear Cahn-Hilliard type models describing diffusion, a thermodynamical framework for the coupling with mechanics is presented. For the mechanical part, besides small deformation theory two different finite deformation elasticity formulations based on elastic Green strain and logarithmic elastic strain, respectively, are introduced and compared. First, a phasefield model for the cathode material *NaxFePO*<sup>4</sup> of NIBs is studied for the first time to understand phase changes, mechanical deformation, and stress evolution in *NaxFePO*<sup>4</sup> electrode particles. As a major novelty, the key parameters in the phase-field model for *NaxFePO*<sup>4</sup> are determined. In this way, our model captures the important feature that distinguishes *NaxFePO*<sup>4</sup> from *LixFePO*4. Furthermore, we study the particle size and average concentration dependent miscibility gap of the nanoscale insertion materials *LixMn*2*O*4, *LixFePO*4, and *NaxFePO*4. Finally, we introduce the nonlocal species concentration theory in terms of diffusion and phase changes from a nonlocal free energy density, compare this theory to the Cahn-Hilliard theory, and show how the nonlocality influences the results.

# **Kurzfassung**

Für viele mobile technische Systeme wie Kommunikation und Elektromobilität werden elektrochemische Energiespeicher benötigt. Lithium-Ionen-Batterien (LIBs) haben in den letzten Jahrzehnten intensive Aufmerksamkeit für die elektrochemische Energiespeicherung auf sich gezogen. Die enorme Nachfrage an LIBs in Kombination mit der begrenzten und ungleichmäßig verteilten Verfügbarkeit von Lithium wird den Lithiumpreis in die Höhe treiben. Im Gegensatz zu LIBs werden Natriumionenbatterien (NIBs), die auf der breiten und reichlichen Verfügbarkeit sowie den niedrigen Kosten von Natrium beruhen, zu einer vielversprechenden Alternative. In den meisten Speichermaterialien treten Phasenänderungen während des regulären Betriebs auf und können nicht vermieden werden. Die jeweiligen Phasen eines Speichermaterials besitzen unterschiedliche Gitterkonstanten, die zu einer Dehnungsfehleranpassung von bis zu einigen hundert Prozent führen. Dadurch entstehen mechanische Spannungen, die zu einer Schädigung der Elektrodenteilchen führt. Vor diesem Hintergrund ist es das Ziel dieser Arbeit, Phasenänderungen und mechanische Spannungen in Elektrodenteilchen mittels Phasenfeldsimulationen zu untersuchen. Thermodynamische Phasenfeldmodelle basieren auf kontinuierlichen Ordnungsparametern, was zu diffusen Grenzflächen zwischen benachbarten Phasen führt, ohne dass die Position mit einer scharfen Grenzfläche aufwändig verfolgt werden muss.

In dieser Arbeit wird für nichtlineare Cahn-Hilliard-Modelle, die die Diffusion beschreiben, ein thermodynamisches Gerüst für die Kopplung mit der Mechanik vorgestellt. Für den mechanischen Teil werden neben der Elastizitätstheorie der kleinen Verzerrungen auch zwei verschiedene Formulierungen der Elastizität für große Verzerrungen, basierend auf dem Green'schen sowie dem logarithmischen Verzerrungsmaß, verglichen. Zum ersten Mal wird ein Phasenfeldmodell für das Kathodenmaterial *NaxFePO*<sup>4</sup> von NIBs untersucht um Phasenänderungen, mechanische Deformation und Spannungsentwicklung in *NaxFePO*<sup>4</sup> -Elektrodenteilchen zu verstehen. Als große Neuheit werden die Schlüsselparameter im Phasenfeldmodell für *NaxFePO*<sup>4</sup> bestimmt. Auf diese Weise erfasst dieses Modell die wichtige Unterscheidung zwischen *NaxFePO*<sup>4</sup> und *LixFePO*4. Darüber hinaus wird die Partikelgröße und die durchschnittliche konzentrationsabhängige Mischungslücke der nanoskaligen Insertionsmaterialien *LixMn*2*O*4, *LixFePO*<sup>4</sup> und *NaxFePO*<sup>4</sup> untersucht. Abschließend wird die nichtlokale Spezieskonzentrationstheorie in Bezug auf Diffusion und Phasenänderungen aus einer nichtlokalen freien Energiedichte eingeführt, mit der Cahn-Hilliard-Theorie verglichen und gezeigt wie die Nichtlokalität die Ergebnisse beeinflusst.

# **Acknowledgment**

The present thesis was written at the Institute for Applied Materials - Materials and Biomechanics (IAM-WBM) of the Karlsruhe Institute of Technology (KIT). I deeply appreciate the financial support from China Scholarship Council (CSC). I would like to express my gratitude to all the people who supported me during my PhD, many of who are not mentioned here.

First I owe my heartiest gratitude to my supervisor Prof. Dr. Marc Kamlah for giving me the possibility to join his research group in this stimulating environment. Especially I am thankful for his inspiring and fruitful discussions, patient explanations, and encouragement during my PhD. This work would never be possible without his support.

My sincere thanks go to the members of my thesis comittee, Prof. Dr. Britta Nestler and Prof. Dr. Robert Stieglitz. I am especially grateful to Prof. Dr. Britta Nestler for accepting to be the second referee, reviewing my thesis and providing valuable suggestions. I would also like to gratefully acknowledge Prof. Dr. Robert Stieglitz for being the chair of the committee.

Furthermore, I would like to thank all my collegues in the IAM-WBM for supporting me whenever I needed help. I am sincerely grateful to Dr. Reiner Mönig, Dr. Dominik Kramer, and Mr. Manfred Janzen for very instructive discussions about the cathode material *NaxFePO*4. I would like to thank Dr. Di Chen for very helpful discussions about crystal structure analysis, as well as the nice support in the whole aspects at the beginning of my PhD. I would like to thank Mr. Oleg Birkholz for his support and valuable suggestions, especially helping me when I had problems by running Linux. I am really grateful to Mr. Friedemann Streich for translating the abstract of this thesis into German, as well as insightful suggestions in the group meeting. I would also appreciate Prof. Dr. Willy Dörfler, Dr. Karthikeyan Chockalingam, and Mr. Fabian Castelli from Institute for Applied and Numerical Mathematics for insightful discussions in the numerical methods. I am especially grateful to Dr. Karthikeyan Chockalingam for extending my knowledge of finite element methods.

In addition, I would like to express my gratitude to all my colleagues and friends for the pleasuring time spent at the IAM-WBM. In particular, I would like to mention Mr. Ewald Ernst, Dr. Kuo Zhang, Dr. Jin Xu, Dr. Dandan Qu, Mr. Yuran Kang, Dr. Zhaoping Luo, Dr. Marigrazia Moscardini, Mr. Likang Luan, Ms. Jin Wang, Ms. Verena Becker, Mr. Junyu Gao, Dr. Lei Chen, Mr. Felix Sutter, Meenakshi, Mr. Huihuang Xia, Dr. Kai Wu, and Dr. Xi Li. Finally, I would like to dedicate a final thank to my parents for their endless love and support. This work could not have been achieved without your love.

# **List of Symbols**

#### Roman letters



#### Greek letters



# **Contents**





# **1 Introduction**

## **1.1 Generalities on battery systems**

Energy conversion and storage have become a growing global concern over the past decades due to urgent energy demand. The renewable energy sources, e.g., solar radiation, wind, and waves, are becoming increasingly popular, although they can be available in intermittent manner [1]. To mitigate this adversity, the chemical energy may help. The chemical energy is the most appropriate form of energy storage in terms of energy density, and batteries provide stored chemical energy with the ability to deliver it as electrical energy with high conversion efficiency and no gaseous exhaust [2].

Depending on the rechargeability, batteries are classified into two categories: the primary and secondary batteries. Primary batteries evolve an irreversible process, converting the chemical energy into electrical energy only once. On the contrary, secondary batteries, also known as rechargeable batteries, are possibly used several times for repeated discharges and charges. This is due to the reversible electrochemical reactions in which the batteries can also be operated in the opposite direction to charge them. Therefore the electrical energy is converted back to chemical energy during charging so that the batteries may be discharged again. We will discuss the reversible process in Section 1.2. However, it should be noticed that even for secondary batteries, capacity fade of batteries occurs after the repeated charging and discharging cycles owing to unavoidable irreversible processes like side reactions, bulk and interfacial resistances of both electrodes and electrolyte, and inhomogenous electrode deformations due to the speices insertion and extraction [3]. Hence, increasing the number of battery cycles is nowadays crucial in order to improve the life time of batteries.

Lithium-ion batteries (LIBs) are widely used and common battery type. As illustrated in Fig. 1.1, LIBs exhibit superb performances in terms of energy density and design flexibility among the existing battery technologies. Indeed, lithium is the lightest metallic element (0.53*gcm*−<sup>3</sup> ) and has a very low redox potential (*E* 0 (*Li*+/*Li*) = −3.04*V* v.s. standard hydrogen electrode), which enables batteries with high voltage and high energy density [1].

Figure 1.1: Comparison of various battery technologies in terms of volumetric and gravimetric energy density [4]. Image reused with the permission of Springer Nature.

Since the commercial LIBs were first announced by Sony in 1991, LIBs have contributed considerably to applications that require portability, such as entertainment, computing and telecommunication equipment. Furthermore, driven by the improved desire for green technologies, the use of batteries has expanded from portable electronics to large scale applications, for example, electric vehicles [5], which would reduce the pollution and secure energy independence.

## **1.2 Design and working principle of lithium-ion batteries**

A battery is composed of several electrochemical cells, connected in series and/or in parallel, to achieve the required voltage and capacity [6]. As depicted in Fig. 1.2, a LIB cell consists of current collectors, anode, cathode, and separator between them. A separator that is commonly made of polymers, in an ideal case, is an ionic conductor and electronic insulator. The ionic conduction of the separator is provided by its pores which are filled with the liquid electrolyte [7]. A current collector, made of copper or aluminum, contacts with the electrode and conducts the electrons to the external circuit.

Figure 1.2: (Left) Schematic of a LIB cell, consisting of current collectors, cathode, anode, and separator. (Right) Illustration of the lithium reduction at the surface of the active particles. The lithium-ions are transported through the electrolyte, which fills the pores in the separator and electrodes. Electronic conduction occurs in the conductive network, made of carbon black and active particles [8].

The electrodes are usually made of a composite of different materials, as shown in Fig. 1.2. The reacting material is formed as small active particles. These mono- or polycrystalline active particles effectively store the lithium-ions [9, 10], and are embedded in the soft and porous binders. A binder, usually made of polyvinylidene difluoride (PVDF), allows to transfer electrolyte through its pores and provides mechanical integrity [11]. To maintain a high electronic conductivity of the electrode, the binder is enriched by so called carbon black, which is a carbon additive of nanosize [7]. An electrically conductive network is built by active particles and carbon additive through which the active particles are connected to the current collector. Insertion materials are usually chosen as the active material. The structure of this material does not change during lithium insertion. As a result, this material acts as a host to provide empty lattice sites that can be occupied by the lithium. For the cathode materials, like lithium manganese oxide *LixMn*2*O*<sup>4</sup> (LMO), lithium cobalt oxide *LixCoO*2, lithium nickel oxide *LixNiO*2, and lithium iron phosphate *LixFePO*<sup>4</sup> (LiFPO), all of them possess their own merits and demerits [9, 10, 12–15]. It should be noticed that LiFPO does not belong to the category of metal oxides. LiFPO has a lower electrochemical potential than the metal oxides, but is thermally more stable [10]. Carbon-based anodes are still very popular anode material in the market today. Alloy anodes such as silicon, aluminum, antimony, and tin are very promising candidates for the next generation of LIBs due to their high theoretical energy density, safe operation potentials, and relatively low cost [7, 16].

The fundamental working principle of a LIB cell is the discrepancy in the electrochemical potential of lithium in two electrodes [17]. Lithium diffuses from the electrode with a high electrochemical potential of lithium to the electrode with a low electrochemical potential of lithium. When the battery is discharged, the difference in the electrochemical potential of lithium of the two electrodes drives lithium-ion to diffuse out of the anode, through the electrolyte, and into the cathode. To keep the electrodes electrically neutral, electrons flow through an external circuit from the anode to the cathode. Both the ionic and the electronic processes are reversed when the battery is charged by an external power source.

## **1.3 Motivation and goals**

The interplay between electrochemistry and mechanics in the battery system has attracted growing scientific interest in the recent years. During charge or discharge, the active material suffers from inhomogeneous volume change, and as a consequence, stresses are induced in the active particles of the electrode. The stresses may lead to fracture of electrode particles. For example, in LIB systems a volume expansion of more than 300% of silicon particles has been observed [18], which can cause particle fracture. Examples for the observation of fracture in storage particles for different materials are shown in Fig. 1.3. Such mechanical degradation leads to capacity loss of a battery over charge and discharge cycles [3]. On the other hand, for thermodynamical reasons, there is a contribution of the stresses to the driving force for diffusion in the host material [19, 20].

Intercalation electrode materials commonly exhibit phase changes during insertion of the intercalation species, such as crystalline silicon (Si) [21], tin, antimony and their oxides [18, 22] for anodes, as well as cathode materials LMO [23] and LiFPO [24–26]. Phase changes may induce large concentration gradients at a mesoscopic scale and thus also large stress magnitudes. Even at low *C*-rates, large stresses may also be induced in the phase-separating electrode materials [27, 28]. For the modeling of phase changes, it is very common to use the Cahn-Hilliard theory where the order parameter is a conserved quantity in order to avoid the need for the cumbersome tracking of the position of a sharp interface [29]. We will review the derivation of the Cahn-Hilliard theory from different approaches in chapter 4.

Figure 1.3: (a) SiSn thin film is cracked after the first few cycles of lithiation and delithiation [30]. Image reused with the permission of Electrochemical Society, Inc. (b) A LiFPO particle is cracked into parts after 60 cycles [31]. Image reused with the permission of Elsevier.

LIBs have attracted intensive attention for electrochemical energy storage over the past decades. The massive use of LIBs combined with the limited and unevenly distributed lithium source will dramatically increase the prices of lithium, and high cost remains a critical problem for the development of LIBs [2]. In contrast to lithium, the wide availability, abundance and low cost of sodium on earth make sodium-ion batteries (NIBs) suitable for large scale energy storage devices in which high energy density becomes less critical [32]. For comparison, properties of lithium and sodium are summarized in Table 1.1. Recently, NIBs have been considered as a promising alternative to LIBs, since sodium and lithium exhibit similar chemical properties so that sodium chemistry could be applied to a similar battery system, and the fundamental principles of the NIBs and LIBs are identical [32]. Just like LIBs, the electrochemical processes in the electrodes of NIBs are also coupled to mechanical properties.


Table 1.1: Sodium versus Lithium characteristics.

Among all cathode materials for NIBs, phosphate based cathode materials are among the best candidates, due to their thermal stability and higher voltage [2]. Olivine *NaxFePO*<sup>4</sup> (NaFPO) has the highest theoretical specific capacity (154*mAhg*−<sup>1</sup> ) compared to the other phosphate based materials (*NaV PO*4*F*, *Na*3*V*2(*PO*4)2*F*<sup>3</sup> and *Na*2*FePO*4*F* etc.) [2], which makes it a promising cathode material for NIBs. Similar to olivine LiFPO, olivine NaFPO also shows phase changes during sodium insertion or extraction [32, 36–39].

As far as we know, no work has been published for the phase-field modeling of the cathode materials of NIBs by now. In this work, a phase-field model for NaFPO of NIBs will be studied for the first time. As a major novelty, the material parameters for NaFPO, in particular, the key parameters in the phase-field model, are determined. In this way, our model captures the important feature that distinguishes NaFPO from LiFPO. Furthermore, the volume expansion of *FePO*<sup>4</sup> upon sodiation to *NaFePO*<sup>4</sup> reaches about 17%, which is much larger than that for LiFPO (about 6.8%) changing from *FePO*<sup>4</sup> to *LiFePO*<sup>4</sup>

[37]. Thus, it is also interesting to study whether or not the small deformation theory (SDT) still has a sufficient capacity to represent the deformation of this material.

On the other hand, the reduction of the electrode particle size down to the nanoscale range can improve the power density [40, 41] and the rate capability [42, 43]. It is found that the miscibility gap between the species-poor phase and the species-rich phase shrinks as the particle size decreases, corresponding to an increase in the mutual solid solubility, and, as a consequence, the misfit strain is reduced [44–51]. The thermodynamics of phase segregation in nanoscale particles is distinctly different from bulk materials, and the particle size dependent miscibility gap plays a nontrivial role in the performance of such nanoscale insertion materials. The previous research works about the particle size miscibility gap have been almost exclusively focused on LiFPO. By now, the particle size and average concentration dependent miscibility gap of the nanoscale insertion materials LMO and NaFPO are still unclear. The question is still open how the mechanical stress affects the particle size and average concentration dependent miscibility gap in intercalation electrode materials. In this work, we will investigate and compare the particle size and average concentration dependent miscibility gap for the three cathode materials LMO, LiFPO, and NaPPO during insertion, using a coupled phase-field model based on the Cahn-Hilliard theory and finite deformation elasticity. We will physically explain the average concentration dependent miscibility gap, and determine the critical particle size below which phase segregation is inhibited for the three cathode materials. We will also study how the mechanical stress affects the particle size and average concentration dependent miscibility gap.

According to Cahn and Hilliard [29], the Cahn-Hilliard theory is derived by truncating the Taylor series of the system free energy, ignoring fourth-order and higher derivatives of the order parameter in the system free energy. Therefore, the system free energy in the Cahn-Hilliard theory depends on the value of the order parameter at a certain position and its immediate neighborhood. As a result, the Cahn-Hilliard theory can not be seen as arising from an exact macroscopic description of microscopic models of interacting particles [52]. In addition, the Cahn-Hilliard equation involves a fourth-order spatial derivative of the order parameter, leading to computational difficulty. In this work, we will introduce a phase-field theory in terms of diffusion and phase changes from a nonlocal free energy density, which can be applied to electrode materials. This theory is named nonlocal species concentration theory (NSC theory). As a special feature, we discuss the nonlocal nature of NSC theory, i.e. in which respect NSC theory accounts for effects in the whole problem domain. Furthermore, the results from this theory and the Cahn-Hilliard theory are compared.

### **1.4 Overview**

The goal of this thesis is to understand phase changes, mechanical deformation, and stress evolution in NaFPO electrode particles of NIBs, and to study the particle size and average concentration dependent miscibility gap for the three cathode materials LMO, LiFPO, and NaPPO. Also, we introduce the NSC theory from a nonlocal free energy density, and show how the nonlocality influences the results.

The work is organized as follows. A literature review about related concepts, methods, and previous studies is given in chapter 2. In chapter 3, fundamentals of continuum mechanics are briefly reviewed. In chapter 4, the derivation of the Cahn-Hilliard theory from various approaches is summarized. Furthermore, the NSC theory is introduced from a nonlocal free energy density. Finally, we discuss the relationship between the CH theory and our theory. In chapter 5, the coupled phase-field model of the Cahn-Hilliard theory and finite deformation elasticity is derived. For the mechanical part, besides SDT two different finite deformation elasticity formulations based on elastic Green strain and logarithmic elastic strain, respectively, are introduced and compared. In chapter 6, a phase-field model for the cathode material NaFPO of NIBs is studied for the first time, and the relevant material parameters of NaFPO are determined. We present and discuss phase changes and stresses in spherical NaFPO electrode particles, and compare the cathode material NaFPO of NIBs to LiFPO of LIBs in terms of phase changes and stresses. In chapter 7, we investigate the particle size and average concentration dependent miscibility gap of three cathode materials LMO, LiFPO, and NaPPO. In chapter 8, we present the results for spherical LMO electrode particles based on the NSC theory, discuss the interface evolution and the nonlocal effect, and compare the results from the NSC theory and the Cahn-Hilliard theory. We conclude our work in chapter 9.

# **2 Literature review**

## **2.1 Phase-field method**

The phase-field method is a powerful computational approach for modeling microstructure evolution in materials at the mesoscale without explicitly tracking the positions of interfaces. It was initiated in 1894 by van der Waals to model a liquid-gas system using a density function that varies continuously at the liquid-gas interface [53]. Approximately sixty years ago, Ginzburg and Landau [54] formulated a phase-field model for superconductivity using a complex valued order parameter and its gradients, and later Cahn and Hilliard [29] carried out a thermodynamic formulation considering the gradients in thermodynamic properties in heterogeneous systems with diffuse interfaces.

Microstructures of physical systems are usually compositional and structural heterogeneous by nature, and can evolve with time. The driving force for microstructural evolution is the possibility to reduce the free energy of the system. The microstructural evolution is described by means of the field variables called order parameters that are continuous functions of time and spatial coordinates. The field variables have different values in different phases to distinguish different microstructures. There are two types of field variables: conserved and nonconserved. The former one has to satisfy the local conservation condition [55, 56]. The temporal and spatial evolutions of the conserved and nonconserved field variables are governed by the well-known Cahn-Hilliard nonlinear diffusion equation as in the present case and the Allen-Cahn relaxation equation, respectively. Both equations can be physically derived from a so-called local microforce balance of Gurtin [57]. The field variables can be an existing physical order parameters such as concentration [29], polarization [58], and magnetization [59]. On the other hand, for the sole purpose of avoiding tracking the interfaces and identifying different phases, dummy fields are introduced as field variables. For example, in the cases of solidification [60] or crack propagation [61], the phase-field that equals zero represents the liquid phase or the unbroken phase while the phase-field with value of unity represents the solid phase or the crack phase.

Phase-field model is based on a diffuse-interface description. The interfaces between phases are continuous across the interfacial regions, as shown in Fig. 2.1a, which is different from the sharp-interface model of the conventional approach in Fig. 2.1b. As a result, an impressive advantage of the phase-field method is that there is no need to follow explicitly the position of the interfaces through mathematical equations during microstructural evolution.

Figure 2.1: (a) Diffuse interface: continuously between phases. (b) Sharp interface: discontinuous at the interface [55]. Image reused with the permission of Elsevier.

The application of the phase-field method has been extensively developed in the past decades, besides solidification [62, 63] and solid-state phase transformations [56, 64], the phase-field method also has been employed in grain growth and coarsening [65–67], the spinodal decomposition of binary mixture [29, 68], crack propagation [61, 69–74], electromigration [75, 76], vesicle membranes in biological applications [77, 78], and multicomponent interdiffusion [79].

## **2.2 Overview of phase-field modeling of electrode materials**

Although multi-physical models have extensively investigated in the past decades, as shown in comprehensive overviews [80–84], it remains a big challenge to simulate the complicated microstructure of electrode materials. The interactions among mechanical, electrochemical, electrical, and thermal fields in the battery system induce complex non linear partial differential equations for the unknown variables (i.e. displacement, concentration, electric potential, and temperature). As a result, it needs more advanced numerical methods to solve such complex multi-physical models. As shown in Figs. 2.2 and 2.3, different processes that occur in a battery during normal operation are described by their corresponding mathematical models. In this work, we will focus on diffusion, mechanical deformation, and phase segregation in active particles at the mesoscale, corresponding to processes of number 3, 4, and 5 in Figs. 2.2 and 2.3, respectively. In what follows, the overview of the phase-field modeling of electrode particles will be presented.

Figure 2.2: Different processes that occur in a battery during normal operation [84]. Image reused with the permission of Springer Nature.

Figure 2.3: Corresponding models for the processes that occur in a battery during normal operation [84]. Image reused with the permission of Springer Nature.

#### **2.2.1 Phase-field modeling of phase changes**

By now, most of the research works about phase-field modeling of phase changes have been focused on LiFPO of LIBs. In 2004, Han et al. [85] used the classical Cahn-Hilliard phase-field model to describe phase segregation in LiFPO electrode particles, and investigated to what extent non-Fickian behavior can affect results from experimental techniques for measuring diffusion coefficients, such as Galvanostatic Intermittent Titration Technique (GITT) and Potentiostatic Intermittent Titration Technique (PITT). According to their results, GITT and PITT, based on the Fick's equation, are still capable of accurately determining diffusion coefficients even if the solid is characterized by a gradient energy term.

In 2008, Singh et al. [86] developed a general continuum theory for phasetransformation dynamics in LiFPO electrode particles by coupling an anisotropic Cahn-Hilliard phase-field model with surface reactions governing the flux of ions across the electrode/electrolyte interface. The general model can induce different transport and phase separation dynamics regimes. For the regime of isotropic bulk-transport-limited phase transformation dynamics, bulk diffusion in all directions is much slower than surface reactions, and the phase boundary is fully contained within the material and moves along the direction of the ionic flux. As shown in Fig. 2.4 a, the shrinking core structure forms. For the regime of anisotropic bulk-transport-limited phase transformation dynamics, as depicted in Fig. 2.4 b, diffusion in the x and z directions is negligible, while diffusion in the more accessible y direction is much slower than surface reactions. The phase boundary still shows a shrinking core in the bulk, but ions are confined to 1D channels in the y direction. Thus, anisotropy alters the shape of the phase interface rather than its basic diffusive dynamics. For the new regime of anisotropic surface-reaction-limited dynamics, surface reactions are much slower than diffusion in the y direction but much faster than diffusion in the x and z directions. As illustrated in Fig. 2.4 c, the phase boundary extends from surface to surface along planes of fast ionic diffusion, rather than forming a classical shrinking core. It should be noticed that these three reaction-diffusion arguments are oversimplified and only give a sense of possible dynamical regimes in the general model.

Figure 2.4: Transport models obtained in different limits of the characteristic timescales for bulk diffusion and surface reactions. Figures show xy cross sections of spherical and platelike single crystals during Li insertion, after phase nucleation has occurred. Lithiated portions of the crystal are shaded, and points outside particles represent flux of Li ions across the electrode/electrolyte interface (shown only for spherical particles). The *FePO*4/*LiFePO*<sup>4</sup> phase boundary is denoted by the dashed line, and arrows indicate movement of the boundary as insertion proceeds. (a) Isotropic bulk transport limited. (b) Anisotropic bulk transport limited. (c) Anisotropic surface reaction limited [86]. Image reused with the permission of Elsevier.

Subsequently, Bai et al. [87] developed a Cahn-Hilliard phase-field model of reaction-limited intercalation in anisotropic LiFPO nanoparticles and predicted that phase separation is suppressed above a critical current. As shown in Fig. 2.5, spinodal decomposition or nucleation leads to moving phase boundaries for small currents. Above a critical current density, the spinodal decomposition is found to disappear, and the particles start to fill homogeneously, which may explain the superior rate capability and long cycle life of nano-LiFPO cathodes.

Figure 2.5: Voltage responses at different constant currents [87]. Image reused with the permission of American Chemical Society.

Recently, Santoki et al. [88] used the classical Cahn-Hilliard equation to study the mesoscopic effect of the surface curvature of the cathodic particle made of LMO. They found that, near the convex region of the particle surface, more sites are available for the applied flux than for the sites that host lithium ions. As a result, the onset of phase segregation prefers to occur in high curvature regions of the particle. Also, the elliptical particle with a higher aspect ratio is subjected to the onset of the phase segregation, prior to the lower ones.

## **2.2.2 Phase-field modeling of diffusion and mechanics**

Here, a related overview of phase-field modeling of diffusion and mechanics is given (see also Zhang and Kamlah [89]). A general phase-field theory for the coupling of the Cahn-Hilliard equation and finite deformation elasticity based on the local dissipation inequality and the so-called local microforce balance has been presented by Gurtin [57]. In this theory, the external microforce itself remains a quantity free to choose, and by taking it equal to zero, the classical Cahn-Hilliard equation is obtained. Walk et al. [28] have identified the influence of the external microforce in the chemical potential for spherical LMO particles. Cogswell and Bazant [50] investigated the effects of elastic coherency strain on the thermodynamics, kinetics, and morphology of intercalation in LiFPO nanoparticles based on a reaction-limited phase-field model. The key parameters in the phase-field model for LiFPO are determined. Their model quantitatively captures the influence of misfit strain on solubility seen in experiments. They concluded that the elastic coherency strain strongly suppresses phase separation during discharge, enhancing the rate capability and extending the cycle life. Tang et al. [90] used a Cahn-Hilliard phase-field model wih SDT to investigate phase separation and crystalline-to-amorphous transformations in spherical isotropic LiFPO nanoparticles, and assessed the conditions under which amorphous phase transitions may occur in LiFPO particles. Huttin and Kamlah [27] considered the coupling between the Cahn-Hilliard equation and SDT for spherical particles of LMO, and demonstrated that large stresses may also occur even at low C-rates.

Anand [91] derived a theory for lithium insertion modeled by the Cahn-Hilliard equation combined with finite deformation elasto-plasticity based on the multiplicative decomposition of the deformation gradient, in which the external microforce is identified as related to the Mandel stress in the chemical potential. To this end, the so-called macroforce and microforce balances are evaluated by following the virtual-power approach of Germain [92] and Gurtin [93]. Subsequently, Di Leo et al. [94] formulated a continuum model which coupled the Cahn-Hilliard-type phase-field theory with finite deformation elasticity based on logarithmic elastic strain. However, Anand's model [91, 94] ignored the fact that logarithmic elastic strain is dependent on the species concentration in the stress chemical potential, although logarithmic elastic strain is also a constitutive variable. As a result, the extremely complicated and important term induced by the concentration dependence of the logarithmic elastic strain in the stress chemical potential is not accounted for. The concentration dependence of the strain tensor in the chemical potential is also confirmed by Gurtin [57]. Interestingly, the external microforce term related to the Mandel stress in Anand's model [91, 94] is just equivalent to the ignored term induced by the concentration dependence of the logarithmic elastic strain, which also is verified by Refs. [28] and [95]. Walk et al. [28] reported that for SDT the identification of the external microforce in Anand's model [91, 94] leads to a doubling of the mechanical coupling term in the chemical potential. Zhao et al. [95] derived a Cahn-Hilliard phase-field model coupled to mechanics based on the Neo-Hookean elasticity without introduction of the external microforce in the chemical potential, and showed that their model agrees with Anand's model [91, 94] which ignored the concentration dependence of the strain tensor in the chemical potential. As remarked by Gurtin [57], no constitutive relation is specified for the external force. Rather, specifying the external microforce by invoking the so-called principle of virtual power may be considered an option. However, in view of the above discussion we do not find it appropriate to adopt this in our work.

## **2.2.3 Studies of the miscibility gap in nanoscale insertion materials**

Here, a related overview of the studies of the miscibility gap in nanoscale insertion materials is given (see also Zhang and Kamlah [96]). The miscibility gap is the difference between the lower and upper bounds of the concentration range of phase segregated states. The research works about the particle size dependent miscibility gap have been almost exclusively focused on LiFPO of LIBs. Meethong et al. [45] observed that the miscibility gap of LiFPO shrinks as the particle size decreases, and they suggested an estimate of the critical particle size (15 *nm*) below which a complete solid solution is achievable at room temperature for LiFPO. Wagemaker et al. [49] revealed that the miscibility gap for LiFPO in small particles not only shrinks, but also strongly depends on the average concentration, which is also called "state of charge" (SOC). In contrast to our common thermodynamic knowledge that the solubility limits are independent of the average concentration, their combined neutron and X-ray diffraction investigation reveals an decreasing miscibility gap that appears to be strongly dependent on the average concentration below particle sizes of 35 *nm*. However, they explained the average concentration dependent solubility limits based on the so-called average solubility limit (i.e. the average compositions in each phase). In our understanding, the average concentration dependent miscibility gap should be related to the minimization of system free energy. What is more, the influence of the elastic strain energy is ignored in their diffuse interface model. The elastic strain energy can suppress phase segregation, and cause larger solid-solution-composition-ranges [27, 50]. Cogswell and Bazant [50] used a coupled phase-field model based on SDT to study the particle size dependent solubility for LiFPO. They performed the phase-field simulations by allowing a square particle at SOC= 0.5 to relax to equilibrium at zero current. Not looking at other SOC values, the average concentration dependent miscibility gap is not accounted for. Welland et al. [51] developed a comprehensive phase-field model based on SDT accounting for facet dependent surface wetting and investigated equilibrium states of LiFPO particles of different size and average lithium concentration, and found that the miscibility gap vanishes for particles of a radius around 5 *nm*, and that the solubility limits change with overall particle lithiation in small particles. However, no flux is applied at the surface of the particles in their work, rather a phase segregated concentration profile or a fluctuation field was used as an initial condition for different surface wetting cases. They have not yet addressed the dynamic loading, which is related to the experimentally more relevant condition of a constant applied flux. Zhao et al. [97] developed a phase-field model accounting for lithium transport in particles, phase separation, and interface reactions across the particle network to study the influence of particle size variation on the performance of LIBs with *V*2*O*<sup>5</sup> nanowires. Their work reveals that both particle size and size variation in electrodes should be small to avoid intra- and inter-particle phase separation.

#### **2.2.4 Nonlocality of phase-field models**

Here, a related overview of nonlocality of phase-field models is given (see also Zhang and Kamlah [98]). Starting with Walter [99], nonlocal models have been studied in the papers [100–111]. In the above mentioned works [99– 107, 109–111], a basic form of the nonlocal energy is used to describe the behavior of materials that exhibit a morphology of phase structures. Compared to the second-order derivative of the order parameter in the system free energy in the Cahn-Hilliard theory, the nonlocal energy form in the nonlocal model has no derivative of the order parameter so that both, discontinuites and oscillations, are allowed. The nonlocal energy is written in an integral form. To be specific, it is related to a weighted integral average of the squared deviation between the order parameter at the point of consideration and the order parameter anywhere else in the problem domain. This will be introduced in Section 4.2.1. In this way, as the system free energy contains an additional term representing an energy average over the entire problem domain, the concept of nonlocality is introduced. The nonlocality is based on the consideration of nonlocal (long-range) particle-particle interactions due to molecular forces from the viewpoint of both classical and modern statistical mechanics [99]. Fosdick and Mason [112] derived a continuum stored energy that includes a nonlocal consideration based on the work of Brandon [103].

Urbachs et al. [113] derived a nonlocal diffuse interface model for the microstructure evolution of tin-lead solder, in which a nonlocal mass fraction is used in order to introduce the nonlocal effect. However, the integral equation related to the normalized weight function in their work [113] does not satisfy a natural weight function property [103, 112, 114]. We will discuss the natural weight function property in Section 4.2.1. As a result, the definition of the nonlocal mass fraction seems to be not physically sound. Indeed, by their expression of what they call the nonlocal mass fraction the nonlocal species concentration, in our wording is not equal to the species concentration for homogenous states, if their theory satisfies the natural weight function property. Also, the nonlocal mass fraction should be explained as a weighted average value of the local mass fraction in the entire problem domain rather than, as they say, over a finite zone in the vicinity of the interface. What is more, the definition of the interface tension coefficient seems to be mathematically inconsistent with the other definitions. Di Leo et al. [94] formulated the gradient micromorphic concentration theory using the principle of virtual power, in which besides conventional species concentration an additional micromorphic concentration is employed. However, this theory is not derived from the nonlocal approach mentioned above, and the relationship between their theory and the Cahn-Hilliard theory seems to not to be clear. Therefore, they just pursued the goal how the solutions from their theory can approach the solutions from the Cahn-Hilliard theory. In particular, they do not discuss aspects of the nonlocal nature of their theory which can account for effects beyond an infinitesimal neighborhood of the point of consideration.

# **3 Basics of continuum mechanics**

### **3.1 Kinematics**

For the description of deformation, a body B is considered in the Euclidean space. A reference configuration is arbitrarily chosen at a fixed reference time. Usually, an undeformed body is preferred to be the reference configuration, and this also corresponds to the initial configuration at time *t* = *t*0. ~*X* denotes an arbitrary material point of B in the reference configuration described by Lagrangian or material coordinates, and~*x* is the spatial position of this material point in the current configuration described by Euler or spatial coordinates. A smooth one-to-one mapping from the material points ~*X* to the spatial point s ~*x* at time *t* describes the motion of B

$$
\vec{\chi} = \vec{\mathcal{X}}(\vec{X}, t). \tag{3.1}
$$

The displacement field

$$
\vec{\mu} = \vec{x} - \vec{X} \tag{3.2}
$$

represents the connection vector between the current and reference configurations. The displacement gradient H is expressed by

$$\mathbf{H} = \mathbf{Grad}\,\vec{\mu}.\tag{3.3}$$

The differentiation of equation (3.1) leads to the deformation gradient

$$\mathbf{F} = \mathbf{Grad}\,\vec{\mathcal{X}} = \mathbf{H} + \mathbf{I},\tag{3.4}$$

which is very crucial in nonlinear continuum mechanics. Here Grad is the gradient operator calculated with respect to the reference configuration, and I is the second order unit tensor. The deformation gradient F is a two-point tensor involving points in two different configurations, and it is a primary measure of deformation [115]. As showin in Fig. 3.1, the fundamental relation between the material line element *d*~*X* and the spatial line element *d*~*x*

$$d\vec{x} = \mathbf{F} \, d\vec{X} \, \tag{3.5}$$

is determined by the deformation gradient.

Figure 3.1: Transformation of a material line element from the reference configuration to the current configuration.

The determinant of the deformation gradient *J* gives the mapping from the infinitesimal reference volume *dV<sup>R</sup>* to the current state *dV*

$$dV = JdV\_R = \det \mathbf{F} \, dV\_R. \tag{3.6}$$

By means of equation (3.6) the following relation is obtained

$$dV = d\vec{A} \cdot d\vec{x} = J d\vec{A}\_R \cdot d\vec{X},\tag{3.7}$$

with *d*~*A* = *dA* ~*n* and *d*~*A<sup>R</sup>* = *dA<sup>R</sup>* ~*n<sup>R</sup>* denoting the infinitesimal surfaces in the current and reference configurations, respectively. Here, ~*n* and ~*n<sup>R</sup>* are the corresponding the unit normal of the surfaces.

Using Equations (3.5) and (3.7) the mapping from the infinitesimal reference area *d*~*A<sup>R</sup>* to the current area *d*~*A* can be obtained as

$$d\vec{A} = J\mathbf{F}^{-\top}d\vec{A}\_R. \tag{3.8}$$

#### **3.2 Stress tensors**

For infinitesimal surfaces,

$$d\mathbf{f} = \mathbf{p}\,dA = \mathbf{P}\,dA\_R \tag{3.9}$$

holds for the infinitesimal force *d*f, where the tractions p and P are defined in the current and reference configurations, respectively, read as

$$\mathbf{p} = \mathbf{T} \ \vec{n} \ \text{and} \ \mathbf{P} = \mathbf{T}\_R \ \vec{n}\_R. \tag{3.10}$$

Here, T*<sup>R</sup>* is the second-order first Piola-Kirchhoff stress tensor, which characterizes the force per unit undeformed area acting on surfaces in the current configuration, and T denotes the symmetric Cauchy or true stress tensor. Combining Equations (3.9) and (3.10) yields

$$\mathbf{T} \, d\vec{A} = \mathbf{T}\_R \, d\vec{A}\_R. \tag{3.11}$$

Using Equation (3.8), the relation between the two stress tensors can be rewritten as

$$\mathbf{T}\_R = J \mathbf{T} \mathbf{F}^{-\top}.\tag{3.12}$$

As shown in Equation (3.12), T*<sup>R</sup>* is not symmetric. Using a pullback operator F −1 , the symmetric second Piola-Kirchhoff stress tensor T˜ *<sup>e</sup>* is introduced

$$\mathbf{\tilde{T}}^{\varepsilon} = \mathbf{F}^{-1} \mathbf{T}\_R = J \mathbf{F}^{-1} \mathbf{T} \mathbf{F}^{-\top},\tag{3.13}$$

and T˜ *<sup>e</sup>* describes the force per unit undeformed area acting on surfaces in the reference configuration.

In small deformation framework, the stretch is small, the rotation is close to unity, and the difference between material and spatial are neglected, i.e the displacement gradient H satisfies

$$|\mathbf{H}| \ll 1,\tag{3.14}$$

and | · | denotes the norm of a vector. The linear strain tensor is defined as the symmetric part of the displacement gradient

$$\boldsymbol{\varepsilon} = \frac{1}{2} (\mathbf{H} + \mathbf{H}^T). \tag{3.15}$$

#### **3.3 Mechanical equilibrium**

An arbitrary body B with the boundary ∂B is subjected to the body force ~*b*, and the traction vector is distributed over its surface. The body is in mechanical equilibrium when the resultant force is zero. The following equations describe the equilibrium conditions in the spatial and material presentations, respectively,

$$\int\_{\partial \mathcal{A}} \mathbf{T} \cdot \vec{n} \, dA + \int\_{\partial \mathcal{B}} \vec{b} \, dV = 0,\tag{3.16}$$

$$
\int\_{\partial \mathcal{A}\_{\mathcal{R}}} \mathbf{T}\_R \cdot \vec{n}\_R \, dA\_R + \int\_{\partial \mathcal{B}\_{\mathcal{R}}} \vec{b} \, dV\_R = 0. \tag{3.17}
$$

Using the Gaussian integral theorem follows

$$\int\_{\mathcal{\partial}} (\operatorname{div} \mathbf{T} + \vec{b}) \, dV = 0,\tag{3.18}$$

$$\int\_{\mathcal{R}\_{\mathcal{R}}} \left( \text{Div} \, \mathbf{T}\_R + \vec{b} \right) \, dV\_R = 0. \tag{3.19}$$

Since the considered body is arbitrary, the above relations are always supposed to hold. Hence, the specified local forms of the above mechanical equilibrium conditions can be deduced

$$\text{div}\,\mathbf{T} + \vec{b} = \mathbf{0},\tag{3.20}$$

$$\text{Div}\,\mathbf{T}\_R + \vec{b} = \mathbf{0}.\tag{3.21}$$

# **4 Phase-field theory of phase changes**

### **4.1 Cahn-Hilliard theory**

The Cahn-Hilliard theory can be physically derived from the thermodynamics by two main approaches. One is based on a variational formulation, following Cahn and Hilliard [29, 68]. The other is derived from the so-called local microforce balance [57]. Both approaches will be reviewed in the following.

#### **4.1.1 System free energy**

To start with, the system free energy is described (see also Zhang and Kamlah [89]). We introduce as an order parameter, depending continuously on space, the species concentration *c*, which is measured in mol per unit volume. The free energy density consisting of two parts is given by

$$
\Psi(\bar{c}, \vec{\nablac}) = \Psi^{m\nu p}(\bar{c}) + \Psi^{\mathcal{S}d}(\vec{\nablac}), \tag{4.1}
$$

where ¯*c* is the dimensionless concentration, which is normalized with respect to the maximum species concentration *cmax* as ¯*c* = *c*/*cmax*, and ~∇ denotes the Nabla operator. The homogeneous free energy density ψ *mwp* is a multiwell potential defining the respective phases. Based on the references [27, 85], the homogeneous free energy density exhibits a double-well structure. The existence of this zone of concavity indicates that homogeneous species concentration states do not always ensure the system free energy to be minimal. In the concavity zone, the system becomes unstable towards phase segregation. The Maxwell construction, which connects the neighborhoods of the two minima by a common tangent, predicts that the system splits into a two phase system and the chemical potentials of both the species-poor phase and species-rich phase are equal. In this sense, the homogeneous free energy density is a multiwell potential where the wells define the respective phases. The second term on the right hand side of Equation (4.1) represents the gradient energy leading to a diffuse interface between two adjacent phases, which is given by

$$\Psi^{gd}(\vec{\nabla c}) = k\_B T\_{ref} N\_A c\_{max} (\frac{1}{2} \lambda \left| \vec{\nabla c} \right|^2). \tag{4.2}$$

Here *k<sup>B</sup>* is the Boltzmann constant, *N<sup>A</sup>* is the Avogadro constant, and *Tre f* is the reference temperature. Furthermore, λ is the gradient energy coefficient with units of length squared controlling the thickness of the diffuse interface. The system free energy of an arbitrary body B of volume *V* is

$$\Psi(\vec{c}) = \int\_{\beta^{\mathcal{J}}} (\Psi^{m\nu p}(\vec{c}) + \Psi^{\mathcal{S}d}(\vec{\nabla}\vec{c}))dV. \tag{4.3}$$

The dimensionless free energy density is introduced as ψ¯ = ψ/(*kBTre fNAcmax*), and any parts of ψ may be normalized in the same way. According to Cahn and Hilliard [29], the system free energy Ψ(*c*¯) can be interpreted physically such that, to the first approximation, the free energy of a small volume of nonuniform solution can be expressed as the sum of two terms, one being the free energy that this volume would have in a homogeneous solution and the other a "gradient energy".

The homogeneous free energy density for a two phase material can be expressed as [27, 85]

$$\begin{aligned} \Psi^{m\nu p}(\bar{c}) &= \quad k\_B T\_{ref} N\_A c\_{max} \left( \alpha\_1 \bar{c} + \frac{\alpha\_2}{2} \bar{c}^2 \right. \\ &\quad + \frac{T}{T\_{ref}} \left( \bar{c} \ln \bar{c} + (1 - \bar{c}) \ln \left( 1 - \bar{c} \right) \right) . \end{aligned} \tag{4.4}$$

The first two terms on the right hand side of Equation (4.4) represent the interaction energy, where positive values of α<sup>1</sup> characterize the energy of inserting a species into the host material, and negative values of α<sup>2</sup> indicate the interaction of neighboring species to be attractive. Concerning α1, it is hard to relate a physical meaning to non-positive values of α1. On the other hand, it has to be noted that α<sup>1</sup> does not occur in the Cahn-Hilliard evolution equation for concentration. According to Huttin [116], there is a critical temperature, namely *T<sup>c</sup>* = −1/4α2*Tre f* , above which the homogeneous free energy density is of convex shape, representing an ideal solution. For temperatures *T* below the critical temperature, the free energy density exhibits a zone of concavity where the homogeneous concentration states are not stable states of the system, and phase segregation becomes possible, as shown in Fig. 4.1. The zone of concavity corresponds to the condition that the second order derivative of the homogeneous free energy density with respect to the concentration is negative. This inequality is never fulfilled if α<sup>2</sup> is not negative. Thus, for a system of noninteracting species, i.e. α<sup>2</sup> = 0, or for a system of species that repel each other, i.e. α<sup>2</sup> > 0, the homogeneous free energy density keeps a convex shape for all *T* > 0. For an attractively interacting species system, i.e. α<sup>2</sup> < 0, depending on the temperature, the homogeneous free energy density may exhibit a zone of concavity. Therefore, at *T* = *Tre f* , the attraction is strong enough to initiate phase segregation for α<sup>2</sup> < −4. The terms multiplied by absolute temperature *T* represent the entropy of mixing, where ¯*c* ln ¯*c* represents the entropy of mixing valid at low concentration, and (1−*c*¯)ln(1−*c*¯) describes the non-dilute solution behavior and represents the entropy of mixing responsible for the saturation effect [95, 117].

Figure 4.1: Normalized homogeneous free energy density versus normalized concentration for different values of absolute temperature *T*.

#### **4.1.2 Derivation from a variational formulation**

The standard derivation of Cahn-Hilliard theory starts from the system free energy, as defined in Equation (4.3), and the variation δΨ with respect to the order parameter *c* reads

$$
\begin{split}
\delta\Psi &= \int\_{\beta\mathcal{B}} (\frac{\partial\Psi^{m\nu p}}{\partial c} \,\delta c + \frac{k\_B T\_{ref} N\_A}{c\_{max}} \lambda \vec{\nabla} c \cdot \delta \vec{\nabla} c) \,dV \\ &= \int\_{\beta\mathcal{B}} (\frac{\partial\Psi^{m\nu p}}{\partial c} - \frac{k\_B T\_{ref} N\_A}{c\_{max}} \lambda \nabla^2 c) \,\delta c \,dV \\ &+ \int\_{\partial\beta\mathcal{B}} \frac{k\_B T\_{ref} N\_A}{c\_{max}} \lambda \vec{\nabla} c \cdot \vec{n} \,\delta c \,dA.\end{split}
\tag{4.5}
$$

Here ∇ 2 is the Laplacian operator. Equation (4.5) yields the natural boundary condition on the entire boundary ∂B

$$
\vec{\nabla}c \cdot \vec{n} = 0.\tag{4.6}
$$

When the interface between phases meets the particle surface, Equation (4.6) enforces that it is perpendicular to the surface [95]. Thus, the chemical potential is defined as the variational derivative of the total free energy with respect to the species concentration

$$
\mu = \frac{\delta \Psi}{\delta c} = \frac{\partial \Psi^{mnp}}{\partial c} - \frac{k\_B T\_{ref} N\_A}{c\_{max}} \lambda \nabla^2 c,\tag{4.7}
$$

and it represents a thermodynamic driving force for diffusion and phase changes. Thus, Equation (4.5) can be rewritten in this form

$$
\delta\Psi = \int\_{\beta\mathcal{\mathcal{S}}} (\frac{\partial \Psi^{\prime \ast \nu \nu}}{\partial c} - \frac{k\_B T\_{ref} N\_A}{c\_{\text{max}}} \lambda \nabla^2 c) \,\delta c \,dV
$$

$$
= \int\_{\beta\mathcal{\mathcal{S}}} \mu \,\,\delta c \,dV.\tag{4.8}
$$

Based on the conservation of mass, the evolution of species transport is given by

$$
\dot{c} = -\operatorname{div}(\vec{J})\tag{4.9}
$$

with the mass flux *J*~ related to the chemical potential µ through the Onsager relation

$$
\vec{J} = -\mathbf{M} \cdot \vec{\nabla}\mu,\tag{4.10}
$$

where the mobility tensor M is non-negative definite. we choose an isotropic mobility according to

$$\mathbf{M}(c) = \mathbf{M}(c)\mathbf{I} \tag{4.11}$$

with the function

$$M(c) = \frac{D\_0 c \left(c\_{\max} - c\right)}{k\_B T\_{ref} N\_A c\_{\max}},\tag{4.12}$$

which is symmetric in the range between zero and maximum concentration. Here, *D*<sup>0</sup> is the diffusion constant.

Combining Equations (4.7), (4.9), (4.10), and (4.11) yields the classical Cahn-Hilliard equation

$$\dot{c} = \text{div}\left(\mathcal{M}(c)\vec{\nabla}\left(\frac{\partial\Psi^{m\nu p}}{\partial c} - \frac{k\_B T\_{ref} N\_A}{c\_{max}}\lambda\nabla^2 c\right)\right). \tag{4.13}$$

#### **4.1.3 Derivation from a microforce balance**

According to Gurtin [57], although the above derivation of Cahn-Hilliard equation is simple and physically sound, it should not be regarded as basic, but rather as precursor of more complete theories. Indeed, such derivation limits the manner in which rate terms enter the equations, and requires a-priori specification of the constitutive equation (4.3). Here, we will review the work of Gurtin [57] to show that how to derive the Cahn-Hilliard equation from a microforce balance.

In the Cahn-Hilliard theory the kinematical process is related to the order parameter *c*, and it is reasonable that there should be "microforces" whose working (expenditure of power) accompanies changes in *c*. Given an arbitrary control volume *R* (subregion of a body B), the system of forces is presumed to be consistent with the microforce balance [57]

$$
\int\_{R} (\pi + \eta) \, dV + \int\_{\partial R} \vec{\tilde{\xi}} \cdot \vec{n} \, dA = 0,\tag{4.14}
$$

where ξ is a vector stress, and π and γ represent, respectively, internal and external forces distributed over the volume of a body B.

The above global law (4.14) should be satisfied for all time and all control volumes *R*, and the local microforce balance is obtained

$$d\text{div}\,\vec{\xi} + \pi + \gamma = 0.\tag{4.15}$$

The general Cahn-Hilliard equation can be derived based on balance of mass, the local microforce balance, and a generalization of the dissipation inequality accommodating diffusion. In order to be consistent with the balance of mass (4.9), here we ignore the external mass supply. Actually, the external mass supply term will be automatically dropped out in the following local dissipation inequality even if it is taken into account.

According to the second law of thermodynamics, the rate at which the free energy of *R* increases can not be more than the working on *R* plus the rate at which free energy is carried into *R* by mass transport. The working is given by

$$\mathcal{W}(\mathsf{R}) = \int\_{\mathsf{R}} \mathsf{y}\dot{\mathsf{c}} \, dV + \int\_{\partial \mathsf{R}} (\vec{\xi} \cdot \vec{n}) \dot{\mathsf{c}} \, dA,\tag{4.16}$$

and the chemical potential through the contributions

$$-\int\_{\partial R} \mu \vec{J} \cdot \vec{n} \, dA \tag{4.17}$$

characterizes the rate at which free energy is carried into *R* by mass transport. Therefore, the appropriate form of the second law of thermodynamics is the dissipation inequality

$$
\int\_{R} \dot{\Psi} \, dV \le \int\_{R} \mathcal{Y} \dot{c} \, dV + \int\_{\partial R} (\vec{\tilde{\xi}} \cdot \vec{n}) \dot{c} \, dA - \int\_{\partial R} \mu \vec{J} \cdot \vec{n} \, dA. \tag{4.18}
$$

The above dissipation inequality can be manipulated further

$$\begin{aligned} &\int\_{R} \dot{\Psi} \, dV - \int\_{R} \gamma \dot{c} \, dV - \int\_{\partial R} (\vec{\tilde{\xi}} \cdot \vec{n}) \dot{c} \, dA + \int\_{\partial R} \mu \vec{J} \cdot \vec{n} \, dA \\ &= \quad \int\_{R} \Psi \, dV - \int\_{R} \gamma \dot{c} \, dV - \int\_{R} \operatorname{div} \left( \dot{c} \vec{\tilde{\xi}} \right) \, dV + \int\_{R} \operatorname{div} \left( \mu \vec{J} \right) \, dV \\ &= \quad \int\_{R} \dot{\Psi} \, dV - \int\_{R} \gamma \dot{c} \, dV - \int\_{R} (\dot{c} \, d\dot{v} \, \vec{\tilde{\xi}} + \vec{\tilde{\xi}} \cdot \vec{\nabla} \dot{c}) \, dV \\ &\quad + \int\_{R} (\mu \, d\dot{v} \, \vec{J} + \vec{J} \cdot \vec{\nabla} \mu) \, dV \\ &\le \quad 0. \end{aligned} \tag{4.19}$$

Using the local microforce balance (4.15) and the balance of mass (4.9), Equation (4.19) yields the local dissipation inequality

$$
\dot{\Psi} - \vec{\xi} \cdot \vec{\nabla} \dot{c} + \vec{J} \cdot \vec{\nabla} \mu + (\pi - \mu) \dot{c} \le 0. \tag{4.20}
$$

According to Gurtin [57], the constitutive equations are considered of the form

$$\begin{aligned} \Psi &= \Psi(c, \vec{\nabla}c, \mu, \vec{\nabla}\mu), & \vec{J} &= \hat{J}(c, \vec{\nabla}c, \mu, \vec{\nabla}\mu), \\ \vec{\xi} &= \hat{\vec{\xi}}(c, \vec{\nabla}c, \mu, \vec{\nabla}\mu), & \pi &= \hbar(c, \vec{\nabla}c, \mu, \vec{\nabla}\mu). \end{aligned} \tag{4.21}$$

All the above constitutive equations should be consistent with the second law of thermodynamics in the form of the local dissipation inequality (4.20) , which leads to the following inequality

$$(\frac{\partial \Psi}{\partial c} + \mathfrak{k} - \mu)\dot{c} + (\frac{\partial \Psi}{\partial \vec{\nabla} c} - \dot{\vec{\xi}}) \cdot \vec{\nabla} \dot{c} + \frac{\partial \Psi}{\partial \mu} \mu + \frac{\partial \Psi}{\partial \vec{\nabla} \mu} \cdot \vec{\nabla} \mu + \vec{J} \cdot \vec{\nabla} \mu \le 0,\tag{4.22}$$

holding for all fields *c* and µ. As a result, we can induce the following relations, which are sufficient for the above inequality (4.22),

$$\begin{aligned} \mathfrak{A}(c, \vec{\nabla}c, \mu) &= \mu - \frac{\partial \Psi(c, \vec{\nabla}c)}{\partial c}, \\ \vec{\xi} &= \hat{\vec{\xi}}(c, \vec{\nabla}c) = \frac{\partial \hat{\Psi}(c, \vec{\nabla}c)}{\partial \vec{\nabla}c}, \end{aligned} \tag{4.23}$$

where ψ and ~ξ are independent of µ and~∇µ, πˆ is independent of ~∇µ, and

$$
\vec{J} \cdot \vec{\nabla} \mu \le 0 \tag{4.24}
$$

for all fields, such that *J*~ has the form of Equation (4.10). Using Equations (4.23) and (4.15) yields

$$
\mu = \frac{\partial \hat{\Psi}(c, \vec{\nabla c})}{\partial c} - d \dot{\nu} (\frac{\partial \hat{\Psi}(c, \vec{\nabla c})}{\partial \vec{\nabla c}}) - \mathcal{y}. \tag{4.25}
$$

It can be interestingly seen that Equation (4.25) with γ = 0 gives the chemical potential obtained from the variation derivative of the total free energy with respect to species concentration, as shown in Equation (4.7).

Substituting Equations (4.10) and(4.25) into the balance of mass (4.9) yields the generalized Cahn-Hilliard equation

$$\dot{c} = \text{div}\left(M(c)\vec{\nabla}\left(\frac{\partial\Psi(c,\vec{\nabla}c)}{\partial c} - d\dot{\nu}(\frac{\partial\Psi(c,\vec{\nabla}c)}{\partial\vec{\nabla}c}) - \gamma\right)\right),\tag{4.26}$$

and the standard Cahn-Hilliard equation is obtained if γ = 0.

### **4.2 Nonlocal species concentration theory**

Here, a nonlocal species concentration theory for diffusion and phase changes is introduced from a nonlocal free energy density (see also Zhang and Kamlah [98]). This theory can be interpreted as an extension of the Cahn-Hilliard theory. In principle, nonlocal effects beyond an infinitesimal neighborhood are taken into account. We first introduce a general formulation of a nonlocal free energy density, and then derive the NSC theory from the nonlocal model.

#### **4.2.1 Nonlocal free energy density**

In the nonlocal model, nonlocal interactions are introduced by the weighted spatial average of the local concentration. Based on [52, 103, 112], the free energy density including truly nonlocal consideration reads as

$$\Psi = \Psi^{m\nu p} \left( \bar{c}(\vec{x}) \right) + \frac{1}{2} A \int\_{\beta \vec{\mathcal{Y}}} \mathcal{o}(\vec{y}; \vec{x}) \left( \bar{c}(\vec{x}) - \bar{c}(\vec{y}) \right)^2 dV. \tag{4.27}$$

Here, ψ *mwp* is the homogeneous free energy density, which has been expressed in Equation (4.4). The second term on the right hand side of Equation (4.27) is the nonlocal free energy density. B (B ⊂ R 3 ) is the entire problem domain containing all material points with position vectors~*x*. In the integral,~*y* denotes the spatial variable for the volume integration over the entire problem domain B. Furthermore, *A* is a material coefficient for the nonlocal particle-particle interaction penalty with units *J*/*m* 3 [112] and ω(~*y*;~*x*) is a positive, symmetric weight function that weights the relative energetic contribution of concentration fluctuations to the free energy density. The nonlocal free energy density in Equation (4.27) vanishes for homogeneous states, i.e. when the concentration at all positions is the same. Thus, the so-called nonlocality becomes active in, say, a phase segregated state, i.e. when the concentration at any position ~*y*, possibly far away from~*x*, is different from the concentration there. The weight function ω(~*y*;~*x*): R <sup>3</sup> → R <sup>+</sup> is a monotonically decreasing function of the distance ρ between ~*x* and ~*y*. It assumes recognizably non-vanishing values over a finite zone in the neighborhood of the material point ~*x*, only, and asymptotically approaches zero as the distance ρ grows to infinity. In this way, the weight function attributes less weight as the distance ρ increases. This satisfies the physical idea that the ability of a material point~*y* to influence the local free energy at a material point~*x* decreases as the distance ρ between the two material points increases. Therefore, the weight function possesses a natural weight function property [103, 112, 114]

$$\int\_{\mathbb{R}^3} \mathbf{o}(\vec{\chi}; \vec{x})dV = 1,\tag{4.28}$$

where, again, ~*y* is the variable of spatial integration. According to Brandon [103], the weight function ω(~*y*;~*x*) ensures that the operator

$$\overline{c} \in L^4(\mathcal{A}; \mathbb{R}) \mapsto \int\_{\mathcal{A}} \mathfrak{o}(\vec{\text{y}}; \vec{\text{x}}) \bar{c}(\vec{\text{y}}) dV \in L^{4/3}(\mathcal{A}; \mathbb{R}) \tag{4.29}$$

is compact, where the space *L* 4 (B;R) is the set of functions ¯*c*

$$
\bar{c}: \mathcal{J} \to \mathbb{R} \tag{4.30}
$$

such that

$$\|\|\tilde{c}\|\|\_{L^{4}(\mathcal{A};\mathbb{R})} \coloneqq \left\{ \int\_{\mathcal{A}^{\sharp}} \tilde{c}^{4}dV \right\}^{1/4} < \infty. \tag{4.31}$$

A similar definition can be applied to the space *L* 4/3 (B;R). For example, the type of the weight function can be chosen of the form [102, 103, 112, 114]

$$\mathfrak{so}(\vec{\mathbf{y}}; \vec{\mathbf{x}}) \sim \frac{e^{-\eta \left| \vec{\mathbf{y}} - \vec{\mathbf{x}} \right|}}{\left| \vec{\mathbf{y}} - \vec{\mathbf{x}} \right|},\tag{4.32}$$

39

which satisfies the natural weight function property and compactness condition when η > 0. The weight function (4.32) is only dependent on the distance ρ = |~*y*−~*x*|. Taking the variation of the free energy (4.27) with respect to *c*(~*x*) to obtain the chemical potential then gives

$$
\mu\_{non} = \frac{\delta\Psi}{\delta c(\vec{\chi})} = \frac{\partial\Psi^{m\nu p}}{\partial c(\vec{\chi})} + \frac{A}{c\_{\text{max}}} \int\_{\beta\theta} a(\vec{\chi};\vec{\chi})(\vec{c}(\vec{\chi}) - \vec{c}(\vec{\chi}))dV,\tag{4.33}
$$

where, similar as in Equation (4.3),

$$
\Psi = \int\_{\beta \overline{\beta}} \Psi(\overline{\mathbf{y}})dV\tag{4.34}
$$

is the system free energy of the entire problem domain B.

We now define two parameters according to

$$
\vec{\mathcal{E}}(\vec{\mathbf{x}}) = \frac{A}{c\_{\text{max}}\mathcal{B}} \int\_{\mathcal{A}} \mathcal{O}(\vec{\mathbf{y}}; \vec{\mathbf{x}}) \vec{c}(\vec{\mathbf{y}}) dV,\tag{4.35}
$$

$$
\beta = \frac{A}{c\_{\text{max}}} \int\_{\beta \overline{\theta}} \mathcal{o}(\vec{y}; \vec{x}) dV. \tag{4.36}
$$

Here ¯*c*ˆ = *c*ˆ/*cmax* is the normalized nonlocal species concentration that is a weighted average value of ¯*c* in the entire problem domain B. The term *A*/(*cmax*β) scales the normalized nonlocal species concentration ¯*c*ˆ such that ¯*c*ˆ is equal to the normalized species concentration ¯*c* for homogenous states. The parameter β is a penalty energy coefficient with units *J*/*mol*, which is related to the penalty energy density, which will be introduced later. It can be normalized as ¯β = β/(*kBTre fNA*).

Using Equations (4.35) and (4.36), the chemical potential can be rewritten to yield

$$
\mu\_{non} = \frac{\delta \Psi}{\delta c(\vec{\chi})} = \frac{\partial \Psi^{m\nu p}}{\partial c(\vec{\chi})} + k\_B T\_{ref} N\_\Lambda \bar{\mathcal{B}} (\bar{c}(\vec{\chi}) - \bar{\hat{c}}(\vec{\chi})).\tag{4.37}
$$

In analogy with Peerlings et al. [114], two different partial differential equations (PDEs) can be derived from the integral formulation (4.35), which will be discussed now in detail.

## **4.2.2 Helmholtz equation governing nonlocal species concentration**

The integral formulation (4.35) can be written in a differential form. First, ¯*c*(~*y*) is expanded into a Taylor series around~*x*

$$\begin{split} \vec{c}(\vec{\uptext{y}}) &= \quad \vec{c}(\vec{\uptext{x}}) + \vec{\nabla}\vec{c}(\vec{\uptext{x}}) \cdot (\vec{\uptext{y}} - \vec{\uptext{x}}) + \frac{1}{2!} \vec{\nabla}\vec{\uptext{V}}\vec{c}(\vec{\uptext{x}}) : (\vec{\uptext{y}} - \vec{\uptext{x}})(\vec{\uptext{y}} - \vec{\uptext{x}}) \\ &+ \frac{1}{3!} \vec{\nabla}\vec{\mathbf{V}}\vec{\uptext{V}}\vec{c}(\vec{\uptext{x}}) \dot{\vec{\uptext{y}}} \langle \vec{\uptext{y}} - \vec{\uptext{x}} \rangle (\vec{\uptext{y}} - \vec{\uptext{x}}) \\ &+ \frac{1}{4!} \vec{\nabla}\vec{\uptext{V}}\vec{\uptext{V}}\vec{\uptext{C}}\vec{c}(\vec{\uptext{x}}) : (\vec{\uptexttext{y}} - \vec{\uptexttext{x}})(\vec{\uptexttext{y}} - \vec{\uptexttext{x}})(\vec{\uptexttexttext{y}} - \vec{\uptexttexttext{x}}) + \cdots \end{split}$$

Again, ~∇ denotes the Nabla operator with respect to spatial position ~*x*. Substituting Equation (4.38) into the integral formulation of the nonlocal species concentration (4.35), yields

¯*c*ˆ(~*x*) = *<sup>A</sup> cmax*β Z B ω(~*y*;~*x*)*c*¯(~*x*)*dV* + *A cmax*β Z B ω(~*y*;~*x*) ~∇*c*¯(~*x*)·(~*y*−~*x*)*dV* + 1 2! *A cmax*β Z B ω(~*y*;~*x*) ~∇~∇*c*¯(~*x*) : (~*y*−~*x*)(~*y*−~*x*)*dV* + 1 3! *A cmax*β Z B ω(~*y*;~*x*) ~∇~∇~∇*c*¯(~*x*) . . .(~*y*−~*x*)(~*y*−~*x*)(~*y*−~*x*)*dV* + 1 4! *A cmax*β Z B ω(~*y*;~*x*) ~∇~∇~∇~∇*c*¯(~*x*) :: (~*y*−~*x*)(~*y*−~*x*)(~*y*−~*x*)(~*y*−~*x*)*dV* +··· . (4.39)

If an isotropic and homogeneous weight function is considered, as the example in Equation (4.32), which is only dependent on the distance ρ = |~*y* −~*x*|, i.e. ω(~*y*;~*x*) = ω(|~*y*−~*x*|) = ω(ρ), odd derivatives vanish upon integration leaving

$$\bar{\hat{c}}(\vec{\mathbf{x}}) = \bar{c}(\vec{\mathbf{x}}) + l^2 \nabla^2 \bar{c}(\vec{\mathbf{x}}) + l\_2^4 \nabla^4 \bar{c}(\vec{\mathbf{x}}) + l\_3^6 \nabla^6 \bar{c}(\vec{\mathbf{x}}) + \dotsb \tag{4.40}$$

with

$$d^2 \quad = \quad \frac{1}{2!} \frac{A}{c\_{\text{max}} \mathfrak{B}} \int\_{\mathcal{A}} \mathfrak{p}^2 \mathfrak{o}(\mathfrak{p}) dV,\tag{4.41}$$

$$d\_2^4 \quad = \quad \frac{1}{4!} \frac{A}{c\_{\text{max}} \mathcal{B}} \int\_{\mathcal{B}} \rho^4 \mathcal{O}(\rho) dV,\tag{4.42}$$

$$I\_3^6 \quad = \underbrace{\frac{1}{6!} \frac{A}{c\_{\max} \mathcal{B}}}\_{\dots \text{-} \atop \cdots \text{-}} \int\_{\mathcal{A}} \rho^6 \alpha(\rho) dV,\tag{4.43}$$

Here the operators ~∇~∇~∇~∇ and ~∇~∇~∇~∇~∇~∇ are simplified to ∇ 4 and ∇ 6 , and so on. Since the series in Equation (4.40) is not truncated, its validity is not limited to some infinitesimal neighborhood of position ~*x*, in principle. In this sense, nonlocal species concentration ¯*c*ˆ is truly nonlocal.

Next, we apply the Laplacian operator to Equation (4.40) and multiply with *l* 2 to obtain

$$d^2\nabla^2\bar{\tilde{c}}(\vec{x}) = l^2\nabla^2\bar{c}(\vec{x}) + l^4\nabla^4\bar{c}(\vec{x}) + l^2l\_2^4\nabla^6\bar{c}(\vec{x}) + \cdots \ . \tag{4.44}$$

Subtracting Equation (4.44) from Equation (4.40) gives

$$\begin{aligned} \tilde{c}(\vec{\boldsymbol{x}}) - l^2 \nabla^2 \tilde{c}(\vec{\boldsymbol{x}}) &= \quad \vec{c}(\vec{\boldsymbol{x}}) + (l\_2^4 - l^4) \nabla^4 \tilde{c}(\vec{\boldsymbol{x}}) \\ &+ (l\_3^6 - l^2 l\_2^4) \nabla^6 \tilde{c}(\vec{\boldsymbol{x}}) + \cdots \, . \end{aligned} \tag{4.45}$$

For vanishing coefficients of the higher order derivatives of ¯*c*, Equation (4.45) reduces to the inhomogenous Helmholtz equation [114]

$$
\bar{\mathcal{E}} - l^2 \nabla^2 \bar{\mathcal{E}} = \bar{\mathcal{C}}.\tag{4.46}
$$

In Appendix A.1, we discuss under which conditions the inhomogenous Helmholtz equation (4.46) becomes the governing equation of nonlocal species concentration ¯*c*ˆ. To this end, we first show that the Green's function of Equation (4.46) satisfies the natural weight function property (4.28). Since additionally the Green's function is isotropic and homogeneous, it is chosen as the weight function for nonlocal species concentration ¯*c*ˆ. Second, we show that ¯*c*ˆ defined in this way satisfies the inhomogenous Helmholtz equation (4.46), indeed.

According to Equation (4.41), the parameter *l* is of the dimension length, and in view of its role in Equation (4.46), we call it characteristic interface length scale. It is taken to be on the order of magnitude of the interface thickness *d* between adjacent phases of the species concentration field ¯*c* [113], and measures the volume in the neighborhood of the material point ~*x* where species concentration ¯*c* contributes significantly to the nonlocal species concentration. For example, for a weight function of the type in Equation (4.32), by Equation (4.41) there is a relation between *l* and the parameter η. Therefore, the characteristic interface length scale *l* is related to the scale of microstructure. Henceforth, we assume *l* to be a constant. The Helmholtz equation (4.46) is a PDE to calculate the nonlocal species concentration. The Helmholtz equation is an elliptic partial differential which represents a time-independent form of the wave equation obtained by separation of variables. It arises, for example, in the study of electromagnetic radiation, seismology, and acoustics.

In order to solve the Helmholtz equation (4.46), we need to introduce a boundary condition. Here the natural boundary condition

$$
\vec{\nabla}\vec{\bar{c}} \cdot \vec{n} = 0\tag{4.47}
$$

is used with~*n* being the unit normal to ∂B. For an interpretation of this boundary condition, integration of Equation (4.46) over the domain B gives

$$
\int\_{\beta\theta} \bar{\mathbf{c}}dV - \int\_{\beta\theta} l^2 \nabla^2 \bar{\mathbf{c}}dV = \int\_{\beta\theta} \bar{\mathbf{c}}dV.\tag{4.48}
$$

The second term of Equation (4.48) vanishes by the divergence theorem and boundary condition (4.47), which gives

$$
\int\_{\partial} \vec{c}dV = \int\_{\partial} \vec{\hat{c}}dV.\tag{4.49}
$$

Thus, the natural boundary condition (4.47) preserves the global species content of the entire problem domain in the nonlocal averaging [114]. In other words, at each time instant, the total content of nonlocal species concentration ¯*c*ˆ is equal to the content of the considered species in the problem domain B.

In summary of the above derivation, together with the boundary condition (4.47), the partial differential equation (4.46) for ¯*c*ˆ is an alternative to the integral representation (4.35) when using the Green's function as weight function. Therefore, the NSC theory based on the Helmholtz equation (4.46) is a special form of the nonlocal model, in which the weight function ω(~*y*;~*x*) is defined as ω(~*y*;~*x*) = *G*(~*y*;~*x*). Indeed, the infinite series of higher-order derivatives of ¯*c* is still implicitly present in the gradient term on the left-hand side of Equation (4.46) [114]. As mentioned before, this indicates that spatial interactions induced by the higher-order derivatives of ¯*c* represent effects at a finite distance and so the NSC theory is truly nonlocal.

### **4.2.3 Decomposition of nonlocal free energy density**

The system free energy

$$\begin{split} \Psi &= \Psi^{m \le p} + \Psi^{penalty} + \Psi^{variance} \\ &= \int\_{\mathcal{B}} \Psi^{m \le p} \left( \bar{\boldsymbol{c}}(\vec{\boldsymbol{y}}) \right) + \Psi^{penalty} \left( \bar{\boldsymbol{c}}(\vec{\boldsymbol{y}}), \bar{\boldsymbol{c}}(\vec{\boldsymbol{y}}) \right) + \Psi^{variance} \left( \bar{\boldsymbol{c}}(\vec{\boldsymbol{y}}) \right) dV \end{split} \tag{4.50}$$

of the NSC theory consists of three parts. For a derivation see Appendix A.2. The first part of Equation (4.50) is the free energy that the entire problem domain would have in a homogenous solution. In the second part, ψ *penalty c*¯(~*x*), ¯*c*ˆ(~*x*) is the penalty energy density of the form

$$\Psi^{penalty}\left(\vec{c}(\vec{x}),\vec{\mathcal{c}}(\vec{x})\right) = \frac{1}{2}c\_{max}\mathcal{B}\left(\vec{c}(\vec{x}) - \vec{\mathcal{C}}(\vec{x})\right)^2. \tag{4.51}$$

This part considers a system penalty energy induced by the the difference between normalized species concentration ¯*c*(~*x*) and normalized nonlocal species concentration ¯*c*ˆ(~*x*). A similar penalty term is also introduced in [94] and [118]. For the last part, we introduce a variance energy density of the form

$$\Psi^{\text{variance}}\left(\vec{\hat{c}}(\vec{\chi})\right) = \frac{1}{2}\mathbf{c}\_{\text{max}}\mathfrak{R}\int\_{\beta\overline{\theta}}\mathfrak{o}(\rho)\left(\vec{\mathfrak{c}}(\vec{\chi}) - \vec{\hat{c}}(\vec{\chi})\right)^{2}dV.\tag{4.52}$$

The last part of Ψ accounts for a system variance energy which depends on the weighted average value of the squared deviation of the local species concentration ¯*c*(~*y*) from normalized nonlocal species concentration ¯*c*ˆ(~*x*). The variance energy density is related to the concept of variance in the probability theory and statistics. It should be noticed that the variance energy density makes no contributuion to diffusion and phase changes due to the independence of the variance energy density with respect to ¯*c*(~*x*). Compared to the free energy density (4.27) in the nonlocal model, the nonlocal free energy density is split into two terms in the NSC theory, i.e. the penalty energy density ψ *penalty* and the variance energy density ψ *variance* .

### **4.3 Comparison of two diffusion theories**

For a better comparison between these two diffusion theories, we first derive the Cahn-Hilliard theory from the nonlocal model. If the terms of order four and higher of Equation (4.40) are neglected, Equation (4.40) yields an approximation of Equation (4.35):

$$
\bar{\hat{c}} = \bar{c} + l^2 \nabla^2 \bar{c}.\tag{4.53}
$$

Using Equation (4.53), Equation (4.37) gives the chemical potential, as shown in Equation (4.7), for the Cahn-Hilliard theory with the gradient energy coefficient λ = ¯β*l* 2 , which is assumed to be a constant, henceforth. According to Equation (4.53), the nonlocal species concentration in a material point depends only on the species concentration and its second-order derivative. As a consequence, the spatial interaction induced by only the second-order derivative of species concentration ¯*c* is limited to the infinitesimal neighborhood of the considered material point. Thus, the Cahn-Hilliard theory is weakly nonlocal.

Compared to the NSC theory, the Cahn-Hilliard theory is a simplified version of the NSC theory, since the NSC theory implicitly accounts for an infinite series of higher-order derivatives of ¯*c* in Equation (4.46) while the Cahn-Hilliard theory only contains the second-order derivative in Equation (4.53). In addition, as will be shown in chapter 8, the interface thickness *d* is controlled by the two parameters ¯β and *l* in the NSC theory, while it depends solely on parameter λ in the Cahn-Hilliard theory. Now we compare the corresponding field equations from the two diffusion theories. Combining Equations (4.9), (4.10), and (4.37) yields the diffusion equation for the NSC theory as

$$\dot{c} = \text{div}\left(M(c)\vec{\nabla}\left(\frac{\partial\Psi^{m\nu p}}{\partial c} + k\_B T\_{ref} N\_A \bar{\beta}\left(\bar{c} - \bar{\hat{c}}\right)\right)\right). \tag{4.54}$$

Here, Equation (4.54) is a partial differential equation (PDE) which involves second-order spatial derivatives in *c*, and the nonlocal species concentration is calculated by solving the Helmholtz equation (4.46). Thus, the NSC theory consists of two coupled second-order PDEs. On the other hand, the Cahn-Hilliard equation (4.13) is a fourth-order PDE. Therefore, the Cahn-Hilliard theory needs stronger continuity requirements on the species concentration than the NSC theory. Instead of a fourth-order equation, the two second-order PDEs of the NSC theory are computationally less demanding.

# **5 Cahn-Hilliard theory with mechanics**

In the previous chapter 4 we introduced the Cahn-Hilliard theory without mechanics. In this chapter, we will derive the coupled model of the Cahn-Hilliard theory and mechanics. We first start with the coupled Cahn-Hilliard theory with SDT. Then, the coupled Cahn-Hilliard theory with finite deformation elasticity is derived, including two finite deformation elasticity formulations based on elastic Green strain and logarithmic elastic strain, respectively. Besides, for the spherically symmetric boundary value problem, the mathematical formulations of the coupled Cahn-Hilliard theories are also derived.

## **5.1 The coupled Cahn-Hilliard theory with small deformation theory**

#### **5.1.1 System free energy**

According to Gurtin [57], in order to take the coupling between diffusion and mechanics into account, the coupling energy or elastic strain energy should be added into the system free energy

$$\begin{split} \Psi(c, \text{grad}c, \mathfrak{e}) &= \int\_{\beta \mathcal{B}} \Psi(c, \text{grad}c, \mathfrak{e}^c)dV \\ &= \int\_{\beta \mathcal{B}} (\Psi^{m \circ p}(c) + \Psi^{\mathfrak{gl}}(\text{grad}c) + \Psi^{cp}(\mathfrak{e}, c))dV. \end{split} \tag{5.1}$$

Here ψ *mwp* and ψ *gd* have been previously referred as the multiwell potential and the gradient energy density, respectively. The coupling energy density ψ *cp* is assumed to be given by

$$\begin{split} \psi^{\varepsilon p} &= \quad \frac{1}{2} \, \mathbf{e}^{\varepsilon} : \mathbb{C} : \mathbf{e}^{\varepsilon} \\ &= \quad G(\mathbf{e}^{\varepsilon} : \mathbf{e}^{\varepsilon} + \frac{\nu}{1 - 2\nu} (\operatorname{tr} \mathbf{e}^{\varepsilon})^2). \end{split} \tag{5.2}$$

Here C is the elasticity tensor

$$\mathbf{C} = G \left( \delta\_{lk} \delta\_{jl} + \delta\_{il} \delta\_{jk} \right) + \frac{2G\nu}{1 - 2\nu} \delta\_{lj} \delta\_{kl},\tag{5.3}$$

which is taken to be constant and isotropic. Accordingly, ν is the Poisson's ratio, and *G* is the shear modulus

$$G = \frac{E}{2(1+\mathbf{v})},\tag{5.4}$$

where *E* is the Young's modulus.

The elastic strain ε *e* is given by

$$
\mathfrak{E}^{\mathfrak{e}} = \mathfrak{E} - \mathfrak{E}^{\mathfrak{s}},
\tag{5.5}
$$

where

$$\mathbf{e}^s = \frac{1}{3}\Omega(c - c\_0)\mathbf{I} \tag{5.6}$$

is the stress-free strain induced by species insertion or extraction, and ε is the total strain tensor that has already introduced in Equation (3.15). Here, Ω is the partial molar volume, and *c*<sup>0</sup> is the initial species concentration.

#### **5.1.2 Linear elasticity**

The Cauchy stress tensor can be derived from the free energy density [57], which leads to the law of linear elasticity

$$\begin{aligned} \mathbf{T} &= \begin{array}{c} \frac{\partial \,\Psi(c, \text{grad}\,c, \mathbf{e})}{\partial \mathbf{e}} \\ = \begin{array}{c} \frac{\partial \,\Psi^{cp}(\mathbf{e}, c)}{\partial \mathbf{e}} \\ = \mathbf{C} : \mathbf{e}^{c}. \end{array} \end{aligned} \tag{5.7}$$

#### **5.1.3 Field equations**

The chemical potential considering the coupling between diffusion and mechanical deformation is a superposition of three terms

$$
\mu = \frac{\delta \Psi}{\delta c} = \mu^{m \circ p} + \mu^{gd} + \mu^{cp} \tag{5.8}
$$

with

$$\begin{aligned} \mu^{\text{mw}p} &= \quad \frac{\partial \,\Psi^{\text{mw}p}}{\partial c} - \underbrace{\text{div}\left(\frac{\partial \,\Psi^{\text{mw}p}}{\partial \,\text{grad}c}\right)}\_{=0} \\ &= \quad k\_B T\_{ref} N\_A \left(\alpha\_1 + \alpha\_2 \bar{c} + \frac{T}{T\_{ref}} \left(\ln \bar{c} - \ln(1 - \bar{c})\right)\right), \tag{5.9} \\ \mu^{\text{gd}} &= \quad \underbrace{\frac{\partial \,\Psi^{\text{gd}}}{\partial c}}\_{=0} - \text{div}\left(\frac{\partial \,\Psi^{\text{gd}}}{\partial \,\text{grad}c}\right) \\ &= \quad -k\_B T\_{ref} N\_A \lambda \,\text{div}\left(\text{grad}c\right), \\ \mu^{\text{cp}} &= \quad \underbrace{\partial \,\Psi^{\text{cp}}}\_{=0} - \underbrace{\text{div}\left(\frac{\partial \,\Psi^{\text{cp}}}{\partial \,\text{grad}c}\right)}\_{=0} \\ &= \quad -\Omega T\_H, \end{aligned} \tag{5.10}$$

where *T<sup>H</sup>* = 1/3*Tii* is the hydrostatic stress. Combining Equations (4.9), (4.10) and (5.8) yields the coupled Cahn-Hilliard diffusion equation

$$\dot{c} = \text{div}\left(\mathcal{M}(c)\vec{\nabla}(\mu^{m\nu p} + \mu^{gd} + \mu^{cp})\right),\tag{5.12}$$

where *M*(*c*) has been defined in Equation (4.12).

It should be noticed that, replacing the term µ *gd* shown in Equation (5.12) by the term

$$
\mu^{penalty} = \frac{\delta \Psi^{penalty}}{\delta c} = k\_B T\_{ref} N\_A \bar{\mathcal{B}}(\bar{c} - \bar{\mathcal{C}}),\tag{5.13}
$$

the coupled NSC diffusion equation with SDT can be obtained.

Finally, based on the balance of momentum, as shown in Equation (3.20), the mechanical equilibrium condition is

$$\text{div}\,\mathbf{T} = \vec{0},\tag{5.14}$$

where the body force is neglected.

Combined with the constitutive equations introduced above, the field equations can be taken as partial differential equations for concentration *c* and displacement vector ~*u*, which need to be solved for given initial and boundary conditions. This is a fourth-order nonlinear initial-boundary-value problem.

### **5.1.4 Mathematical formulation of the spherically symmetric boundary value problem**

In this part, the phase segregation problem of species insertion into or extraction from a cathodic particle is mathematically formulated. In order to avoid costly three-dimensional simulations, here, a spherical particle of radius *R*<sup>0</sup> under spherically symmetric boundary conditions is considered, and the spherical coordinate system (*r*,θ,ϕ) is introduced. It is assumed that not only the particle geometry, but also the unknown variables, such as the species concentration *c* and the displacement vector ~*u*, as well as the boundary conditions holding at the electrode particle surface are invariant under rotation with respect to both the θ and ϕ coordinates [116]. As a result, the three-dimensional problem is allowed to be replaced by an equivalent one-dimensional problem, as sketched in Fig. 5.1. Due to the spherical symmetry, all fields are expressed as a function of the time *t* and the radial coordinate 0 ≤ *r* ≤ *R*0:

$$c = c(r, t),\tag{5.15}$$

$$
\vec{\mu} = \mu\_r(r, t)\vec{e}\_r. \tag{5.16}
$$

Figure 5.1: Spherical particle of the spherically symmetry: the three-dimensional problem reduces to a one-dimensional problem.

In numerical implementations, the scaling of the equations by material parameters may lead to larger numerical errors, especially when the prefactors are very small or large. Therefore, the boundary value problem is implemented in a normalized form. In this work, the relevant normalized parameters are listed in Table 5.1.


Table 5.1: Normalized parameters.

For the spherically symmetric boundary value problem considered here, the strain tensor is diagonal, i.e. entries only exist on the main diagonal. The radial and tangential components of it are given by (referring to Equation (3.15))

$$
\varepsilon\_r = \frac{\partial u\_r}{\partial r},\tag{5.17}
$$

$$
\mathfrak{e}\_{\mathfrak{t}} = \begin{array}{c} \underline{u}\_{\mathfrak{t}} \\ \underline{r} \end{array} \tag{5.18}
$$

According to Equation (5.5), the elastic strain ε *e* is expressed as

$$\begin{array}{ll} \mathbf{c}^{\varepsilon} &=& \mathbf{c} - \mathbf{c}^{\varepsilon} \\ &=& \begin{pmatrix} \frac{\partial \boldsymbol{u}\_{\varepsilon}}{\partial \boldsymbol{r}} - \frac{\boldsymbol{\Omega}}{\boldsymbol{\mathcal{S}}} (\boldsymbol{c} - \boldsymbol{c}\_{0}) & \mathbf{0} & \mathbf{0} \\\\ \mathbf{0} & \frac{\boldsymbol{u}\_{\varepsilon}}{\boldsymbol{r}} - \frac{\boldsymbol{\Omega}}{\boldsymbol{\mathcal{S}}} (\boldsymbol{c} - \boldsymbol{c}\_{0}) & \mathbf{0} \\\\ \mathbf{0} & \mathbf{0} & \frac{\boldsymbol{u}\_{\varepsilon}}{\boldsymbol{r}} - \frac{\boldsymbol{\Omega}}{\boldsymbol{\mathcal{S}}} (\boldsymbol{c} - \boldsymbol{c}\_{0}) \end{pmatrix}\_{\{\vec{\boldsymbol{r}}\_{\varepsilon}, \vec{\boldsymbol{r}}\_{\theta}, \vec{\boldsymbol{r}}\_{\theta}\}}. \end{array} . \tag{5.19}$$

Furthermore, based on the law of linear elasticity (5.7), the components of the stress tensor are given by

$$\begin{aligned} T\_r &= \,^2 2G \left[ \mathfrak{e}\_r + \frac{\nu}{1 - 2\nu} \left( \mathfrak{e}\_r + 2\mathfrak{e}\_t \right) \right] \\ &- \frac{1}{3} \left( 1 + \frac{3\nu}{1 - 2\nu} \right) \mathfrak{A} (c - c\_0) \right], \\\\ T\_t &= \,^2 2G \left[ \mathfrak{e}\_t + \frac{\nu}{1 - 2\nu} \left( \mathfrak{e}\_r + 2\mathfrak{e}\_t \right) \\ &- \frac{1}{3} \left( 1 + \frac{3\nu}{1 - 2\nu} \right) \mathfrak{A} (c - c\_0) \right]. \end{aligned} \tag{5.20}$$

Based on Equation (5.14), the mechanical equilibrium condition of the spherically symmetry is formulated as

$$\frac{\partial T\_r}{\partial r} + \frac{2}{r}(T\_r - T\_l) = 0.\tag{5.22}$$

Now we derive the coupled diffusion equation of the spherically symmetric boundary value problem. For the chemical potential, µ *mwp* has been expressed in Equation (5.9). According to Equations (5.10) and (5.11), µ *gd* and µ *cp* are expressed as

$$\begin{split} \mu^{\otimes d} &= -k\_B T\_{ref} N\_A \lambda \operatorname{div} (\operatorname{grad} \bar{c}) \\ &= -k\_B T\_{ref} N\_A \lambda (\frac{\partial^2 \bar{c}}{\partial r^2} + \frac{2}{r} \frac{\partial \bar{c}}{\partial r}), \\ \mu\_{\text{tot}} &= \operatorname{conv} \end{split} \tag{5.23}$$

$$\begin{aligned} \mu^{cp} &= -\mathfrak{\Omega}T\_H\\ &= -G\left[\frac{2}{3}\mathfrak{\Omega}^2(c-c\_0) - \frac{2}{3}\mathfrak{\Omega}\frac{\partial u\_r}{\partial r} - \frac{4}{3}\mathfrak{\Omega}\frac{u\_r}{r} \\ &- \frac{2\mathfrak{v}\mathfrak{\Omega}}{1-2\mathfrak{v}}\left(-\mathfrak{\Omega}(c-c\_0) + \frac{\partial u\_r}{\partial r} + 2\frac{u\_r}{r}\right)\right]. \end{aligned} \tag{5.24}$$

55

The gradient of the chemical potential is the driving force for diffusion, and the components of it are obtained as:

$$\operatorname{grad} \mu^{\text{mv}p} \quad = \quad \frac{k\_B T\_{ref} N\_A}{\vec{c} (1 - \vec{c})} (1 + \mathfrak{a}\_2 \vec{c} (1 - \vec{c})) \frac{\partial \vec{c}}{\partial r} \vec{e}\_r,\tag{5.25}$$

$$\operatorname{grad} \mu^{cd} \quad = \begin{array}{c} \ -k\_B T\_{ref} N\_A \lambda \left( \frac{\partial^3 \bar{c}}{\partial r^3} + \frac{2}{r} \frac{\partial^2 \bar{c}}{\partial r^2} - \frac{2}{r^2} \frac{\partial \bar{c}}{\partial r} \right) \vec{e}\_r, \\\ . \end{array} \tag{5.26}$$

$$\begin{split} \text{grad}\,\mu^{cp} &= \quad G\left[\frac{2}{3}\Omega^2 \frac{\partial c}{\partial r} - \frac{2}{3}\Omega \frac{\partial^2 u\_r}{\partial r^2} - \frac{4}{3}\Omega \left(\frac{1}{r}\frac{\partial u\_r}{\partial r} - \frac{u\_r}{r^2}\right) \\ &- \frac{2\nu\Omega}{1-2\nu} \left(-\Omega \frac{\partial c}{\partial r} + \frac{\partial^2 u\_r}{\partial r^2} + 2\left(\frac{1}{r}\frac{\partial u\_r}{\partial r} - \frac{u\_r}{r^2}\right)\right)\right] \vec{e}\_r. \end{split} \tag{5.27}$$

Substituting the above Equations (5.25)-(5.27) into Equation (5.12), we can derive the final dimensionless coupled diffusion equation

$$\begin{split} 0 &= \quad \vec{r}^2 \frac{\partial \vec{c}}{\partial \vec{r}} + \frac{\partial}{\partial \vec{r}} \left[ -\vec{r}^2 (1 + \alpha\_2 \vec{c} (1 - \vec{c})) \frac{\partial \vec{c}}{\partial \vec{r}} \\ &+ \vec{r}^2 \frac{\lambda}{R\_0^2} \vec{c} (1 - \vec{c}) \left( \frac{\partial^3 \vec{c}}{\partial \vec{r}^3} + \frac{2}{\vec{r}} \frac{\partial^2 \vec{c}}{\partial \vec{r}^2} - \frac{2}{\vec{r}^2} \frac{\partial \vec{c}}{\partial \vec{r}} \right) \\ &- \vec{r}^2 \vec{c} (1 - \vec{c}) \vec{G} \left[ \frac{2}{3} \vec{\Omega}^2 \frac{\vec{\partial} \vec{c}}{\partial \vec{r}} - \frac{2}{3} \vec{\Omega}^2 \frac{\partial^2 \vec{u}\_r}{\partial \vec{r}^2} - \frac{4}{3} \vec{\Omega} \left( \frac{1}{\vec{r}} \frac{\partial \vec{u}\_r}{\partial \vec{r}} - \frac{\vec{u}\_r}{\vec{r}^2} \right) \right. \\ &- \frac{2 \nu \vec{\Omega}}{1 - 2 \nu} \left( -\vec{\Omega} \frac{\partial \vec{c}}{\partial \vec{r}} + \frac{\partial^2 \vec{u}\_r}{\partial \vec{r}^2} + 2 \left( \frac{1}{\vec{r}} \frac{\partial \vec{u}\_r}{\partial \vec{r}} - \frac{\vec{u}\_r}{\vec{r}^2} \right) \right) \right]. \end{split} (5.28)$$

We can find that the coupled Cahn-Hilliard diffusion equation involves fourthorder spatial derivatives in the concentration and third-order spatial derivatives in the radical displacement. Here, if the term related to shear modulus *G*¯ is not taken into account, the above coupled diffusion equation (5.28) will degenerate to the pure Cahn-Hilliard diffusion equation.

Figure 5.2: Boundary and symmetry conditions for the spherically symmetric boundary value problem.

We now discuss the boundary and initial conditions for the spherically symmetric boundary value problem (see also Zhang and Kamlah [89]). The boundary conditions are sketched in Fig. 5.2. At the surface, the particle is assumed to be stress free,

$$\mathbf{T}(R\_0, t) \cdot \vec{n} = \vec{0},\tag{5.29}$$

where ~*n* refers to the outgoing unit vector normal to the particle surface. We choose a spatially independent mass flux at the surface as

$$\vec{J}(R\_0, t) \cdot \vec{n} = \begin{cases} -\frac{C c\_{\text{max}} R\_0}{3600 \cdot 3} & \text{for } \ c(R\_0, t) \le c\_{\text{max}} \\ 0 & \text{for } \ c(R\_0, t) = c\_{\text{max}} \end{cases} . \tag{5.30}$$

Here *C* is the *C*-rate, and *C* = *n* means that the amount of a species of a fully charged particle would flow into the particle within 1/*n* hours. Once the maximum concentration *cmax* is reached at the surface, the mass flux will be stopped. Here neglecting surface wetting [119], the natural boundary condition is expressed as

$$
\vec{a}\,\text{grad}\,c\cdot\vec{n} = \frac{\partial c}{\partial r}(\mathcal{R}\_0, t) = 0.\tag{5.31}
$$

57

In addition, at the particle center, the boundary conditions

$$
\mu\_r(0, t) = 0,\tag{5.32}
$$

$$\frac{\partial c}{\partial r}(0,t) = 0,\tag{5.33}$$

$$\frac{\partial^3 c}{\partial r^3}(0, t) = 0\tag{5.34}$$

are to be satisfied. Here, Equation (5.32) guarantees the continuity of the particle at the center, and Equations (5.32)-(5.34) are needed to ensure the spherical symmetry. What is more, these three boundary conditions ensure the physical requirement that the flux at the particle center vanishes [116]. Finally, the initial conditions are given by

$$
\mu\_r(r,0) = 0,\tag{5.35}
$$

$$c(r,0) = 0.\tag{5.36}$$

## **5.2 The coupled Cahn-Hilliard theory with finite deformation elasticity**

#### **5.2.1 Kinematics**

The kinematics developed here is related to the basic fields shown in Table 5.2, some of which have been discussed in chapter 3.1 (see also Zhang and Kamlah [89, 96]).


Table 5.2: Basic fields in the kinematics.

Consider a motion~*x* =~χ(~*X*,*t*) of a homogeneous isotropic body B. The kinematical basis for the coupling between finite deformation elasticity and diffusion is the multiplicative decomposition [91]

$$\mathbf{F} = \mathbf{F}^{\epsilon} \mathbf{F}^{s},\tag{5.37}$$

as illustrated in Fig. 5.3. Here F *e* is the elastic deformation gradient that is a mapping from the intermediate configuration onto the current configuration, and the concentration induced deformation gradient F *s* caused by species insertion or extraction is a mapping from the reference configuration onto the intermediate configuration. F *s* is assumed to be purely volumetric in the form

$$\mathbf{F}^s = \beta^s \mathbf{I} = \sqrt[3]{1 + \mathbf{\Omega}(c\_R - c\_0)} \mathbf{I},\tag{5.38}$$

where *c<sup>R</sup>* is the species concentration, measured in mol per unit reference volume.

Using Equations (3.6) and (5.37) yields

$$J = J^{\epsilon} J^{s},\tag{5.39}$$

where *J <sup>e</sup>* = *det*F e and *J <sup>s</sup>* = *det*F <sup>s</sup> = 1+Ω(*c<sup>R</sup>* −*c*0).

Figure 5.3: Multiplicative decomposition of F [89].

The right and left polar decompositions of F *e* are expressed by

$$\mathbf{F}^{\epsilon} = \mathbf{R}^{\epsilon} \mathbf{U}^{\epsilon} = \mathbf{V}^{\epsilon} \mathbf{R}^{\epsilon},\tag{5.40}$$

where R *e* is a rotation, and U *e* and V *e* are the symmetric, positive-definite elastic right and left stretch tensors:

$$\mathbf{U}^{\epsilon} = \sqrt{\mathbf{F}^{\epsilon \top} \mathbf{F}^{\epsilon}} \text{ and } \mathbf{V}^{\epsilon} = \sqrt{\mathbf{F}^{\epsilon} \mathbf{F}^{\epsilon \top}}. \tag{5.41}$$

In addition, the right elastic Cauchy-Green tensor is

$$\mathbf{C}^{\epsilon} = \mathbf{U}^{\epsilon 2} = \mathbf{F}^{\epsilon \top} \mathbf{F}^{\epsilon}. \tag{5.42}$$

The elastic Green strain tensor is then defined as

$$\begin{aligned} \mathbf{E}^{\epsilon} &=& \frac{1}{2} \left( \mathbf{C}^{\epsilon} - \mathbf{I} \right) \\ &=& \frac{1}{2} \left( \mathbf{F}^{\epsilon \top} \mathbf{F}^{\epsilon} - \mathbf{I} \right) \\ &=& \frac{1}{2} \left( \left( 1 + \boldsymbol{\Omega} (c\_{\mathcal{R}} - c\_{0}) \right)^{-\frac{2}{3}} \mathbf{F}^{\top} \mathbf{F} - \mathbf{I} \right) \end{aligned} \tag{5.43}$$

Next, the spectral representations of U *e* and V *e* are

$$\mathbf{U}^{e} = \sum\_{a=1}^{3} \lambda\_a^{e} \vec{r}\_a \otimes \vec{r}\_a \text{ and } \mathbf{V}^{e} = \sum\_{a=1}^{3} \lambda\_a^{e} \left(\mathbf{R}^{e} \vec{r}\_a\right) \otimes \left(\mathbf{R}^{e} \vec{r}\_a\right), \tag{5.44}$$

respectively where (~*r*1,~*r*2,~*r*3) and (λ *e* 1 ,λ *e* 2 ,λ *e* 3 ) are the orthonormal eigenvectors and the positive eigenvalues of U *e* respectively. According to Anand [91], the logarithmic elastic strain tensor and the spatial logarithmic elastic strain tensor can be obtained as

$$\mathbf{E}\_{\log}^{e} = \ln \mathbf{U}^{e} = \sum\_{a=1}^{3} \ln \lambda\_{a}^{e} \vec{r}\_{a} \otimes \vec{r}\_{a} \,\mathrm{,}\tag{5.45}$$

61

$$\begin{array}{rcl} \mathbf{E}\_{\log,H}^{\epsilon} & = & \ln \mathbf{V}^{\epsilon} \\ & = & \sum\_{a=1}^{3} \ln \lambda\_{\alpha}^{\epsilon} \left( \mathbf{R}^{\epsilon} \vec{r}\_{\alpha} \right) \otimes \left( \mathbf{R}^{\epsilon} \vec{r}\_{a} \right) . \end{array} \tag{5.46}$$

Using Equations (5.37), (5.38), (5.41), and (5.45), the above logarithmic elastic strain tensor can be rewritten as

$$\begin{aligned} \mathbf{E}\_{\log}^{\epsilon} &= \quad \ln \mathbf{U}^{\epsilon} \\ &= \quad \ln \sqrt{\mathbf{F}^{\epsilon \top} \mathbf{F}^{\epsilon}} \\ &= \quad \ln(\frac{1}{\sqrt[3]{1 + \Omega(c\_R - c\_0)}} \sqrt{\mathbf{F}^{\top} \mathbf{F}}) \\ &= \quad \ln \sqrt{\mathbf{F}^{\top} \mathbf{F}} - \frac{1}{3} \ln(1 + \Omega(c\_R - c\_0)) \mathbf{I}. \end{aligned} \tag{5.47}$$

Thus, the logarithmic elastic strain tensor involves an eigenstrain accounting for the strain developed by composition changes from the reference state.

#### **5.2.2 System free energy**

The system free energy of the coupled theory, in the reference volume, includes three parts [57, 91]

$$\begin{split} \Psi\_{\mathcal{R}} &= \int\_{\beta \mathcal{B}\_{\mathcal{R}}} \Psi \mathcal{R} \, dV\_{\mathcal{R}} \\ &= \int\_{\beta \mathcal{B}\_{\mathcal{R}}} \left( \Psi\_{\mathcal{R}}^{m \nu p} (c\_{\mathcal{R}}) + \Psi\_{\mathcal{R}}^{gd} (\text{Grad} \, c\_{\mathcal{R}}) + \Psi\_{\mathcal{R}}^{cp} \right) \, dV\_{\mathcal{R}}, \end{split} \tag{5.48}$$

where ψ *mwp R* and ψ *gd R* have been previously discussed as the multiwell potential and the gradient energy density, respectively. Here, entities with subscript *R* indicate these to be defined in the reference configuration. The additional term ψ *cp R* is the coupling energy density defining the coupling between diffusion and mechanics. As shown in Fig. 5.3, the elastic Green strain tensor is defined in the intermediate configuration. Thus, the coupling energy of the entire system in terms of elastic Green strain tensor in the intermediate configuration is given by

$$
\Psi\_{l}^{cp} = \int\_{\beta\mathcal{B}\_{\mathcal{J}}} \frac{1}{2} \,\mathbf{E}^{\varepsilon} : \mathbf{C} : \mathbf{E}^{\varepsilon} dV\_{l}. \tag{5.49}
$$

According to the multiplicative decomposition of F, we can obtain

$$dV\_I = \det \mathbf{F}^s \, dV\_R = J^s \, dV\_R. \tag{5.50}$$

Using Equations (5.49) and (5.50) yields the coupling energy in terms of elastic Green strain tensor in the reference configuration

$$
\Psi\_{\mathcal{R}}^{cp} = \int\_{\beta \mathcal{B}\_{\mathcal{R}}} \frac{1}{2} \, J^s \, \mathbf{E}^e \, :\, \mathbf{C} : \mathbf{E}^e dV\_{\mathcal{R}}.\tag{5.51}
$$

Therefore, the coupling energy density in terms of elastic Green strain tensor in the reference configuration is given by

$$\begin{aligned} \left(\boldsymbol{\Psi}\_{\boldsymbol{\kappa}}^{cp}(\boldsymbol{c}\_{R}, \mathbf{E}^{\epsilon})\right) &=& \frac{1}{2} \, J^{s} \, \mathbf{E}^{\epsilon} : \mathbf{C} : \mathbf{E}^{\epsilon} \\ &=& G J^{s} \left[\mathbf{E}^{\epsilon} : \mathbf{E}^{\epsilon} + \frac{\nu}{1 - 2\nu} (\text{tr}\, \mathbf{E}^{\epsilon})^{2}\right], \end{aligned} \tag{5.52}$$

Similarly, the coupling energy density in terms of logarithmic elastic strain tensor in the reference configuration is derived as

$$\begin{aligned} \left(\Psi\_{\kappa}^{cp}(c\_{\mathbb{R}}, \mathbf{E}\_{\log}^{\epsilon})\right) &=& \frac{1}{2} J^s \, \mathbf{E}\_{\log}^{\epsilon} : \mathbf{C} : \mathbf{E}\_{\log}^{\epsilon} \\ &=& G J^s \left[\mathbf{E}\_{\log}^{\epsilon} : \mathbf{E}\_{\log}^{\epsilon} + \frac{\mathbf{v}}{1 - 2\mathbf{v}} \left(\operatorname{tr} \mathbf{E}\_{\log}^{\epsilon}\right)^2\right]. \end{aligned} \tag{5.53}$$

#### **5.2.3 Elasticity based on elastic Green strain**

According to Gurtin [57], the thermodynamical potential relation for the first Piola-Kirchhoff stress T*<sup>R</sup>* is derived from the local dissipation inequality as

$$\mathbf{T}\_R = \frac{\partial \,\Psi\_R}{\partial \mathbf{F}}.\tag{5.54}$$

The second Piola-Kirchhoff stress T˜ *<sup>e</sup>* related to intermediate configuration is symmetric considering the symmetry of the Cauchy stress T:

$$\mathbf{\tilde{T}}^{e} = J^{e} \mathbf{F}^{e-1} \mathbf{T} \mathbf{F}^{e-\top}. \tag{5.55}$$

Using Equations (3.12), (5.54) and (5.55), we can obtain the relation

$$
\tilde{\mathbf{T}}^{\epsilon} = \frac{2}{J^s} \frac{\partial \Psi\_R}{\partial \mathbf{C}^{\epsilon}} = \frac{1}{J^s} \frac{\partial \Psi\_R}{\partial \mathbf{E}^{\epsilon}}.\tag{5.56}
$$

Using Equations (5.48) and (5.52), Equation (5.56) then yields the isotropic elasticity law based on elastic Green strain as

$$\begin{split} \tilde{\mathbf{T}}^{\varepsilon} &= \quad \frac{1}{J^{s}} \left( \underbrace{\frac{\partial \,\Psi\_{\tilde{R}}^{\mathrm{mvp}}}{\partial \mathbf{E}^{\varepsilon}}}\_{=0} + \underbrace{\frac{\partial \,\Psi\_{\tilde{R}}^{\mathrm{gsd}}(\mathrm{Grad}\,c\_{R})}{\partial \mathbf{E}^{\varepsilon}}}\_{=0} + \underbrace{\frac{\partial \,\Psi\_{\tilde{R}}^{\mathrm{cp}}(\mathrm{Grad}\,c\_{R})}{\partial \mathbf{E}^{\varepsilon}}}\_{=0} \right) \\ &= \quad \frac{1}{J^{s}} \frac{\partial \,\Psi\_{\tilde{R}}^{\mathrm{cp}}(c\_{R}, \mathbf{E}^{\varepsilon})}{\partial \mathbf{E}^{\varepsilon}} \\ &= \quad 2G\mathbf{E}^{\varepsilon} + 2G\frac{\mathbf{V}}{1 - 2\mathbf{V}}(\mathrm{tr}\,\mathbf{E}^{\varepsilon})\mathbf{I} \\ &= \quad 2G\mathbf{E}^{\varepsilon} + \left(K - \frac{2}{3}G\right)(\mathrm{tr}\,\mathbf{E}^{\varepsilon})\mathbf{I}, \end{split} \tag{5.57}$$

where the bulk modulus *K* is expressed by

$$K = \frac{E}{\Im(1 - 2\nu)}.\tag{5.58}$$

#### **5.2.4 Elasticity based on logarithmic elastic strain**

In the following, we will introduce another isotropic elasticity law based on logarithmic elastic strain (see also Zhang and Kamlah [89]). According to Anand [91], the second Piola-Kirchhoff stress T˜ *<sup>e</sup>* is derived from the free energy density according to

$$
\tilde{\mathbf{T}}^{\epsilon} = \frac{2}{J^s} \frac{\partial \Psi\_R}{\partial \mathbf{C}^{\epsilon}}.\tag{5.59}
$$

Here, Equation (5.59) agrees with Equation (5.56) which is derived from the local dissipation inequality of Gurtin [57]. The Mandel stress related to the intermediate configuration is defined by

$$\mathbf{M}^{\epsilon} = \mathbf{C}^{\epsilon} \tilde{\mathbf{T}}^{\epsilon}.\tag{5.60}$$

Using Equations (5.42) and (5.45) yields

$$\mathbf{E}\_{\log}^{\epsilon} = \ln(\mathbf{C}^{\epsilon \frac{1}{2}}) . \tag{5.61}$$

Hence, using Equations (5.59) and (5.61), Equation (5.60) then gives

$$\mathbf{M}^{\epsilon} = \mathbf{C}^{\epsilon} \frac{2}{J^{s}} \frac{\partial \,\Psi\_{\mathsf{R}}}{\partial \mathbf{C}^{\epsilon}} = \frac{1}{J^{s}} \frac{\partial \,\Psi\_{\mathsf{R}}}{\partial \mathbf{E}\_{\log}^{\epsilon}}.\tag{5.62}$$

Using Equations (5.48) and (5.53), Equation (5.62) then yields the isotropic elasticity law based on logarithmic elastic strain as

$$\begin{split} \mathbf{M}^{e} &= \quad \frac{1}{J^{s}} \left( \underbrace{\frac{\partial \Psi\_{\mathrm{R}}^{\mathrm{mw}p}}{\partial \mathbf{E}\_{\log}^{\epsilon}}}\_{=0} + \underbrace{\frac{\partial \Psi\_{\mathrm{R}}^{\mathrm{gd}}(\mathrm{Grad}\,\mathrm{c}\_{R})}{\partial \mathbf{E}\_{\log}^{\epsilon}}}\_{=0} + \frac{\partial \Psi\_{\mathrm{R}}^{\mathrm{e}p}(\mathrm{c}\_{R}, \mathbf{E}\_{\log}^{\epsilon})}{\partial \mathbf{E}\_{\log}^{\epsilon}} \right) \\ &= \quad \frac{1}{J^{s}} \frac{\partial \Psi\_{\mathrm{R}}^{\mathrm{e}p}(\mathrm{c}\_{R}, \mathbf{E}\_{\log}^{\epsilon})}{\partial \mathbf{E}\_{\log}^{\epsilon}} \\ &= \quad \, 2G\mathbf{E}\_{\log}^{\epsilon} + \left( K - \frac{2}{3}G \right) \left( \mathrm{tr}\,\mathbf{E}\_{\log}^{\epsilon} \right) \mathbf{I}. \end{split} \tag{5.63}$$

Furthermore, Equations (5.42), (5.55) and (5.60) yields the corresponding Cauchy stress

$$\mathbf{T}^{\top} = \frac{1}{J^{\epsilon}} \mathbf{R}^{\epsilon} \mathbf{M}^{\epsilon} \mathbf{R}^{\epsilon^{\top}}.\tag{5.64}$$

#### **5.2.5 Field equations**

The chemical potential considering the coupling between diffusion and mechanics is given by

$$
\mu\_{\mathcal{R}}(c\_{\mathcal{R}}, \mathbf{F}) = \frac{\delta \Psi\_{\mathcal{R}}}{\delta c\_{\mathcal{R}}} = \mu\_{\mathcal{R}}^{m \circ p} + \mu\_{\mathcal{R}}^{\otimes d} + \mu\_{\mathcal{R}}^{cp}, \tag{5.65}
$$

where µ *mwp R* and µ *gd R* have been expressed in Equations (5.9) and (5.23), respectively. It should be noticed that, *c* and *r* in Equations (5.9) and (5.23) are replaced by *c<sup>R</sup>* and *R*, respectively, indicating these defined in the reference configuration.

For the coupling chemical potential µ *cp R* , we consider two forms of it for the two different finite deformation elasticities:

$$\begin{split} \mu\_{\kappa}^{cp}(c\_{R}, \mathbf{E}^{\varepsilon}) &= \quad \frac{\partial \Psi\_{R}^{cp}}{\partial c\_{R}} - \underbrace{\mathrm{Div}\left(\frac{\partial \Psi\_{R}^{cp}}{\partial \mathrm{Grad}c\_{R}}\right)}\_{=0} \\ &= \quad \frac{\partial \Psi\_{R}^{cp}}{\partial c\_{R}} \\ &= \quad \frac{\partial \Psi\_{R}^{cp}}{\partial J^{s}} \frac{\partial J^{s}}{\partial c\_{R}} + \frac{\partial \Psi\_{R}^{cp}}{\partial \mathbf{E}^{\varepsilon}} : \frac{\partial \mathbf{E}^{\varepsilon}}{\partial c\_{R}} \\ &= \quad \Omega(\frac{1}{2} \, \mathbf{E}^{\varepsilon} : \mathbf{C} : \mathbf{E}^{\varepsilon}) + J^{s} \mathbf{C} : \mathbf{E}^{\varepsilon} : \frac{\partial \mathbf{E}^{\varepsilon}}{\partial c\_{R}} \\ &= \quad \Omega(\frac{1}{2} \, \mathbf{E}^{\varepsilon} : \mathbf{C} : \mathbf{E}^{\varepsilon}) - \frac{1}{3J^{s}} \Omega \frac{\partial \Psi\_{R}^{cp}}{\partial \mathbf{E}^{\varepsilon}} \mathbf{C}^{\varepsilon} : \mathbf{I} \\ &= \quad \Omega(\frac{1}{2} \, \mathbf{E}^{\varepsilon} : \mathbf{C} : \mathbf{E}^{\varepsilon}) - \frac{1}{3} \Omega \, \mathbf{tr} \mathbf{M}^{\varepsilon}, \end{split} \tag{5.66}$$

µ *cp R* (*cR*,E *e log*) = ∂ψ*cp R* ∂ *c<sup>R</sup>* −Div ∂ψ*cp R* ∂ Grad*c<sup>R</sup>* | {z } =0 = ∂ψ*cp R* ∂ *c<sup>R</sup>* = ∂ψ*cp R* ∂ *J s* ∂ *J s* ∂ *c<sup>R</sup>* + ∂ψ*cp R* ∂E *e log* : ∂E *e log* ∂ *c<sup>R</sup>* = Ω( 1 2 E *e log* : C : E *e log*) +*J <sup>s</sup>*C : E *e log* : ∂E *e log* ∂ *c<sup>R</sup>* = Ω( 1 2 E *e log* : C : E *e log*)− 1 3*J s* Ω ∂ψ*cp R* ∂E *e log* : I = Ω( 1 2 E *e log* : C : E *e log*)− 1 3 ΩtrM*<sup>e</sup>* . (5.67)

Here, the term multiplied by <sup>∂</sup><sup>E</sup> *e* ∂ *cR* or ∂E *e log* ∂ *cR* represents the concentration dependence of the strain tensor in the chemical potential. As shown in the above derivation, this term is equivalent to the external microforce term related to the Mandel stress in Anand's model [91, 94].

Combining Equations (4.9), (4.10), and (5.65) yields the coupled Cahn-Hilliard diffusion equation

$$\dot{c}\_{R} = \text{Div}\left(M(c\_{R})\,\text{Grad}(\mu\_{\text{g}}^{m\nu p} + \mu\_{\text{g}}^{gd} + \mu\_{\text{g}}^{cp})\right),\tag{5.68}$$

where *M*(*c*) has been defined in Equation (4.12). Again, replacing the term µ *gd R* by the term µ *penalty*, the coupled NSC diffusion equation with finite deformation elasticity can be obtained.

Finally, based on the balance of momentum, as shown in Equation (3.21), the mechanical equilibrium condition is

$$\operatorname{div} \mathbf{T}\_{\mathbf{R}} = \vec{0},\tag{5.69}$$

where the body force is neglected.

## **5.2.6 Mathematical formulations of the spherically symmetric boundary value problem**

#### **The coupled diffusion equation with elastic Green strain**

For the spherically symmetric boundary value problem, as mentioned in section 5.1.4, only the radial component of displacement field is non-zero. As a result, the displacement gradient H expressed by Equation (3.3) can be rewritten in a diagonal form, only containing two different components *H*<sup>11</sup> and *H*<sup>22</sup> = *H*<sup>33</sup>

$$\mathbf{H} = \mathbf{Grad}\,\vec{u}(R,t) = \begin{pmatrix} \frac{\partial u\_R}{\partial R} & 0 & 0\\ 0 & \frac{u\_R}{R} & 0\\ 0 & 0 & \frac{u\_R}{R} \end{pmatrix}\_{\{\vec{x}\_R,\vec{x}\_\theta,\vec{x}\_\theta\}}.\tag{5.70}$$

Again, entities with subscript *R* indicate those defined in the reference configuration. According to Equation (3.4), the deformation gradient F is expressed as

$$\mathbf{F} = \mathbf{H} + \mathbf{I} = \begin{pmatrix} \frac{\partial u\_R}{\partial R} + 1 & 0 & 0 \\\\ 0 & \frac{u\_R}{R} + 1 & 0 \\\\ 0 & 0 & \frac{u\_R}{R} + 1 \end{pmatrix}\_{\{\vec{x}\_R, \vec{x}\_\theta, \vec{x}\_\theta\}}.\tag{5.71}$$

Combining Equations (5.37), (5.38), and (5.71) yields the elastic deformation gradient

$$\mathbf{F}^{e} = \frac{1}{\mathcal{B}^{s}} \mathbf{F} = \frac{1}{\mathcal{B}^{s}} \begin{pmatrix} \frac{\partial u\_{R}}{\partial R} + 1 & 0 & 0\\ 0 & \frac{u\_{R}}{R} + 1 & 0\\ 0 & 0 & \frac{u\_{R}}{R} + 1 \end{pmatrix}\_{\{\vec{\mathcal{E}}\boldsymbol{\vec{\mathcal{R}}}\_{\mathcal{E}}, \vec{\mathcal{E}}\boldsymbol{\varrho}\}},\tag{5.72}$$

where β *<sup>s</sup>* = p<sup>3</sup> 1+Ω(*c<sup>R</sup>* −*c*0) (see Equation (5.38)). Accordingly, based on Equations (5.43) and (5.72), the elastic Green strain tensor is expressed as

$$\begin{array}{rcl} \mathbf{E}^{\epsilon} &=& \frac{1}{2} \left( \mathbf{F}^{\epsilon \top} \mathbf{F}^{\epsilon} - \mathbf{I} \right) \\\\ &=& \begin{pmatrix} \frac{(\frac{d\boldsymbol{w}\_{R}}{d\boldsymbol{R}} + 1)^{2}}{2\boldsymbol{\beta}^{\epsilon 2}} - 0.5 & \mathbf{0} & \mathbf{0} \\\\ \mathbf{0} & \frac{(\frac{\mathbf{w}\_{R}}{d\boldsymbol{\beta}^{\epsilon 2}} + 1)^{2}}{2\boldsymbol{\beta}^{\epsilon 2}} - 0.5 & \mathbf{0} \\\\ \mathbf{0} & \mathbf{0} & \frac{(\frac{\mathbf{w}\_{R}}{d\boldsymbol{\beta}^{\epsilon 2}} + 1)^{2}}{2\boldsymbol{\beta}^{\epsilon 2}} - 0.5 \end{pmatrix}\_{\{\vec{x}\_{R}, \vec{x}\_{\theta}, \vec{x}\_{\theta}\}} \\\ & & \tag{5.73} \end{array} \tag{5.73}$$

Substituting Equation (5.73) into Equation (5.52), the coupling energy density in terms of elastic Green strain tensor is given by

$$\begin{split} \Psi\_{\mathbb{R}}^{cp}(c\_{R}, \mathbf{E}^{c}) &= \quad GJ^{s} \left[ \mathbf{E}^{c} : \mathbf{E}^{c} + \frac{\nu}{1 - 2\nu} (\operatorname{tr} \mathbf{E}^{c})^{2} \right] \\ &= \quad G \left( 1 + (c\_{R} - c\_{0}) \Omega \right) \left[ \frac{1}{4} \left( \frac{\left( 1 + \frac{\partial u\_{R}}{\partial R} \right)^{2}}{\left( 1 + (c\_{R} - c\_{0}) \Omega \right)^{\frac{2}{3}}} - 1 \right)^{2} \\ &\quad + \frac{1}{2} \left( \frac{\left( 1 + \frac{u\_{R}}{R} \right)^{2}}{\left( 1 + (c\_{R} - c\_{0}) \Omega \right)^{\frac{2}{3}}} - 1 \right)^{2} \\ &\quad + \frac{\nu}{1 - 2\nu} \left( \frac{1}{2} \left( \frac{\left( 1 + \frac{\partial u\_{R}}{\partial R} \right)^{2}}{\left( 1 + (c\_{R} - c\_{0}) \Omega \right)^{\frac{2}{3}}} - 1 \right) \\ &\quad + \frac{\left( 1 + \frac{u\_{R}}{R} \right)^{2}}{\left( 1 + \left( c\_{R} - c\_{0} \right) \Omega \right)^{\frac{2}{3}}} - 1 \right)^{2} \right]. \end{split}$$

Using Equation (5.55), the relationship between the second Piola-Kirchhoff stress tensor T˜ *<sup>e</sup>* and the Cauchy stress tensor T can be rewritten as

$$\mathbf{T} = \frac{1}{J^{e}} \mathbf{F}^{e} \tilde{\mathbf{T}}^{e} \mathbf{F}^{e^{\top}}.\tag{5.75}$$

Combining Equations (3.12), (5.57), and (5.75), the first Piola-Kirchhoff stress tensor is given by

$$\begin{aligned} \mathbf{^rT\_R} &= \begin{aligned} \frac{J}{J^\epsilon} \mathbf{^r \tilde{T}^\epsilon} \mathbf{F}^\epsilon \mathbf{F}^{\epsilon^\top} \mathbf{F}^{-\top} \\ &= \begin{aligned} \mathcal{B}^\delta \mathbf{F} \tilde{\mathbf{T}}^\epsilon \\ &= \begin{array}{c} 2G\beta^s \mathbf{F} \left( \mathbf{E}^\epsilon + \left( K - \frac{2}{3} G \right) (\text{tr} \, \mathbf{E}^\epsilon) \mathbf{I} \right) . \end{aligned} \end{aligned} \tag{5.76}$$

As shown in Equation (5.76), we see the major diagonal entries of T*<sup>R</sup>* are nonzero, and *TR*<sup>22</sup> = *TR*<sup>33</sup> due to the spherical symmetry. The components of it are therefore given by

$$\begin{aligned} T\_{R\_{11}} &= \quad G \left( 1 + (c\_R - c\_0) \Omega \right)^{\frac{1}{3}} \left( \frac{\partial u\_R}{\partial \mathcal{R}} + 1 \right) \\ &\quad \cdot \left[ \frac{\left( \frac{\partial u\_R}{\partial \mathcal{R}} + 1 \right)^2}{\left( 1 + (c\_R - c\_0) \Omega \right)^{\frac{2}{3}}} - 1 \\ &\quad + \frac{\nu}{1 - 2\nu} \left( \frac{\left( \frac{\partial u\_R}{\partial \mathcal{R}} + 1 \right)^2}{\left( 1 + (c\_R - c\_0) \Omega \right)^{\frac{2}{3}}} + \frac{\left( 2 \frac{u\_R}{R} + 1 \right)^2}{\left( 1 + (c\_R - c\_0) \Omega \right)^{\frac{2}{3}}} - 3 \right) \right], \end{aligned} \tag{5.77}$$

$$\begin{array}{rcl} T\_{R2} &=& G \left( 1 + \left( c\_{\mathcal{R}} - c\_{0} \right) \mathfrak{Q} \right)^{\frac{1}{3}} \left( \frac{u\_{\mathcal{R}}}{\mathcal{R}} + 1 \right) \\\\ & \cdot \left[ \frac{\left( \frac{\partial u\_{\mathcal{R}}}{\partial \mathcal{R}} + 1 \right)^{2}}{\left( 1 + \left( c\_{\mathcal{R}} - c\_{0} \right) \mathfrak{Q} \right)^{\frac{2}{3}}} - 1 \\\\ & + \frac{\mathbf{v}}{1 - 2\mathbf{v}} \left( \frac{\left( \frac{\partial u\_{\mathcal{R}}}{\partial \mathcal{R}} + 1 \right)^{2}}{\left( 1 + \left( c\_{\mathcal{R}} - c\_{0} \right) \mathfrak{Q} \right)^{\frac{2}{3}}} + \frac{\left( 2 \frac{u\_{\mathcal{R}}}{\mathcal{R}} + 1 \right)^{2}}{\left( 1 + \left( c\_{\mathcal{R}} - c\_{0} \right) \mathfrak{Q} \right)^{\frac{2}{3}}} - 3 \right) \end{array} \tag{5.78}$$

Furthermore, due to spherical symmetry, the mechanical equilibrium (5.69) is expressed as

$$\frac{\partial T\_{R\_{11}}}{\partial R} + \frac{2}{R} (T\_{R\_{11}} - T\_{R\_{22}}) = 0. \tag{5.79}$$

µ

Now we derive the coupled diffusion equation of the spherically symmetric boundary value problem. First, the coupled chemical potential µ *cp R* (*cR*,E *e* ), as shown in Equation (5.66), is expressed as

*cp R* (*cR*,E *e* ) = *G*Ω 1 4 <sup>−</sup>1<sup>+</sup> 1+ ∂*uR* ∂*R* 2 (1+ (*c<sup>R</sup>* −*c*0)Ω) 2 3 2 + 1 2 −1+ 1+ *uR R* 2 (1+ (*c<sup>R</sup>* −*c*0)Ω) 2 3 !2 + ν 1−2ν <sup>−</sup>1<sup>+</sup> 1 2 <sup>−</sup>1<sup>+</sup> 1+ ∂*uR* ∂*R* 2 (1+ (*c<sup>R</sup>* −*c*0)Ω) 2 3 + 1+ *uR R* 2 (1+ (*c<sup>R</sup>* −*c*0)Ω) 2 3 !# +*G*(1+ (*c<sup>R</sup>* −*c*0)Ω) −Ω 1+ ∂*uR* ∂*R* 2 3(1+ (*c<sup>R</sup>* −*c*0)Ω) 5 3 · <sup>−</sup>1<sup>+</sup> 1+ ∂*uR* ∂*R* 2 (1+ (*c<sup>R</sup>* −*c*0)Ω) 2 3 − 2Ω 1+ *uR R* 2 3(1+ (*c<sup>R</sup>* −*c*0)Ω) 5 3 −1+ 1+ *uR R* 2 (1+ (*c<sup>R</sup>* −*c*0)Ω) 2 3 !

$$\begin{split} &+\frac{2\nu}{1-2\nu} \\ &\cdot \left(\frac{-\mathfrak{Q}\left(1+\frac{\partial u\_{R}}{\partial R}\right)^{2}}{3\left(1+(c\_{R}-c\_{0})\mathfrak{Q}\right)^{\frac{2}{3}}}-\frac{2\mathfrak{Q}\left(1+\frac{u\_{R}}{R}\right)^{2}}{3\left(1+(c\_{R}-c\_{0})\mathfrak{Q}\right)^{\frac{2}{3}}}\right) \\ &\cdot \left(-1+\frac{1}{2}\left(-1+\frac{\left(1+\frac{\partial u\_{R}}{\partial R}\right)^{2}}{\left(1+(c\_{R}-c\_{0})\mathfrak{Q}\right)^{\frac{2}{3}}}\right) \\ &+\frac{\left(1+\frac{u\_{R}}{R}\right)^{2}}{\left(1+\left(c\_{R}-c\_{0}\right)\mathfrak{Q}\right)^{\frac{2}{3}}}\right). \end{split} \tag{5.80}$$

Substituting Equations (5.9), (5.23), and (5.80) into Equation (5.68), we can derive the dimensionless coupled diffusion equation

$$\begin{split} 0 &= \quad \mathcal{R}^2 \frac{\partial \bar{c}\_R}{\partial \bar{t}} + \frac{\partial}{\partial \mathcal{R}} \left[ \frac{\partial}{\partial \mathcal{R}} \left( -\mathcal{R}^2 \mathcal{M}(c\_{\mathcal{R}}) (\mu\_{\mathcal{R}}^{m \nu p} + \mu\_{\mathcal{R}}^{gd} + \mu\_{\mathcal{R}}^{cp} (c\_{\mathcal{R}}, \mathbf{E}^{\epsilon})) \right) \right]. \end{split} \tag{5.81}$$

For the detailed derivation of the final dimensionless coupled diffusion equation, please see Appendix A.3.1.

The boundary and initial conditions for the spherically symmetric boundary value problem are the same as those of the previous Cahn-Hilliard diffusion with SDT, as sketched in Fig. 5.4. Again, entities with subscript *R* in Fig. 5.4 indicate those defined in the reference configuration.

Figure 5.4: Boundary and symmetry conditions for the spherically symmetric boundary value problem in the reference configuration [89].

#### **The coupled diffusion equation with logarithmic elastic strain**

According to Equations (5.41) and (5.72), the right stretch tensor is expressed as

$$\begin{array}{rcl} \mathbf{U}^{e} &=& \sqrt{\mathbf{F}^{e \top} \mathbf{F}^{e}} \\ &=& \frac{1}{\mathcal{B}^{s}} \mathbf{F} \\ &=& \left(1 + (c\_{R} - c\_{0})\boldsymbol{\Omega}\right)^{-\frac{1}{\frac{1}{\beta}}} \begin{pmatrix} \frac{\partial u\_{R}}{\partial R} + 1 & 0 & 0 \\ 0 & \frac{u\_{R}}{R} + 1 & 0 \\ 0 & 0 & \frac{u\_{R}}{R} + 1 \end{pmatrix}\_{\{\ddot{x}\_{\overline{R}, \ddot{x}\_{\theta}, \ddot{x}\_{\theta}\}} \\ & & & \tag{5.82} \end{array} . \tag{5.82}$$

Therefore, combining the polar decomposition F *<sup>e</sup>* = R *e*U *e* and the above Equation (5.82) yields the rotation tensor

$$\mathbf{R}^e = \mathbf{I}.\tag{5.83}$$

Furthermore, the spectral decomposition of U *e* can be obtained as

$$\begin{aligned} \mathbf{U}^{\epsilon} &= \sum\_{a=1}^{3} \lambda\_a^{\epsilon} \vec{r}\_a \otimes \vec{r}\_a \\ &= \lambda\_1^{\epsilon} \vec{e}\_{\mathbb{R}} \otimes \vec{e}\_{\mathbb{R}} + \lambda\_2^{\epsilon} (\mathbf{I} - \vec{e}\_{\mathbb{R}} \otimes \vec{e}\_{\mathbb{R}}), \end{aligned} \tag{5.84}$$

where

$$\mathcal{A}\_{1}^{\epsilon} \quad = \quad (1 + (c\_R - c\_0)\Omega)^{-\frac{1}{3}} (\frac{\partial u\_R}{\partial R} + 1), \tag{5.85}$$

$$
\lambda\_{\frac{\nu}{2}}^{\epsilon} = \left(1 + (c\_R - c\_0)\Omega\right)^{-\frac{1}{3}} (\frac{\mu\_R}{R} + 1). \tag{5.86}
$$

Thus, the logarithmic elastic strain tensor is expressed as

$$\begin{array}{rclcl} \mathbf{E}\_{\log}^{\epsilon} &=& \ln \mathbf{U}^{\epsilon} \\ &=& \sum\_{\alpha=1}^{3} \ln \lambda\_{\alpha}^{\epsilon} \vec{r}\_{\alpha} \otimes \vec{r}\_{\alpha} \\ &=& \begin{pmatrix} \ln(\frac{1}{\vec{B}^{\epsilon}}(\frac{\partial u\_{\mathsf{R}}}{\partial \vec{R}}+1)) & 0 & 0 \\ 0 & \ln(\frac{1}{\vec{B}^{\epsilon}}(\frac{u\_{\mathsf{R}}}{\vec{R}}+1)) & 0 \\ 0 & 0 & \ln(\frac{1}{\vec{B}^{\epsilon}}(\frac{u\_{\mathsf{R}}}{\vec{R}}+1)) \end{pmatrix}\_{\begin{pmatrix} \vec{x}\_{\mathsf{R}} \vec{x}\_{\mathsf{0}}, \vec{x}\_{\mathsf{0}} \end{pmatrix} \\ & (5.87) \end{array} . \tag{5.87}$$

Here, the spatial logarithmic elastic strain tensor E *e log*,*H* is equal to the logarithmic elastic strain tensor E *e log* due to R *<sup>e</sup>* = I. Substituting Equation (5.87) into Equation (5.53), the coupling energy density in terms of logarithmic elastic strain tensor is given by

$$\begin{split} \boldsymbol{\psi}\_{R}^{cp}(\boldsymbol{c}\_{R},\mathbf{E}\_{log}^{\epsilon}) &= \quad GJ^{s}\left[\mathbf{E}\_{log}^{\epsilon}:\mathbf{E}\_{log}^{\epsilon}+\frac{\mathbf{v}}{1-2\mathbf{v}}\left(\text{tr}\,\mathbf{E}\_{log}^{\epsilon}\right)^{2}\right] \\ &= \quad G\left(1+(c\_{R}-c\_{0})\Omega\right)\left[\left(\ln\left(\frac{\frac{\partial u\_{R}}{\partial R}+1}{\left(1+(c\_{R}-c\_{0})\Omega\right)^{\frac{1}{3}}}\right)\right)^{2}\right] \\ &\quad + 2\left(\ln\left(\frac{\frac{u\_{R}}{R}+1}{\left(1+(c\_{R}-c\_{0})\Omega\right)^{\frac{1}{3}}}\right)\right)^{2} \\ &\quad + \frac{\mathbf{v}}{1-2\mathbf{v}}\left(\ln\left(\frac{\frac{\partial u\_{R}}{\partial R}+1}{\left(1+(c\_{R}-c\_{0})\Omega\right)^{\frac{1}{3}}}\right)\right) \\ &\quad + 2\ln\left(\frac{\frac{u\_{R}}{R}+1}{\left(1+(c\_{R}-c\_{0})\Omega\right)^{\frac{1}{3}}}\right)\bigg)^{2}\bigg]. \end{split} \tag{5.88}$$

Combining Equations (3.12), (5.64), (5.63) and (5.83), the first Piola-Kirchhoff stress tensor is given by

$$\begin{aligned} \mathbf{T}\_R &= \begin{array}{c} \frac{J}{J^\epsilon} \mathbf{M}^\epsilon \mathbf{F}^{-\top} \\\\ &= \begin{array}{c} J^\epsilon \left( 2G \mathbf{E}\_{\log}^{\epsilon} + \left( K - \frac{2}{3} G \right) \left( \text{tr} \, \mathbf{E}\_{\log}^{\epsilon} \right) \mathbf{I} \right) \mathbf{F}^{-\top} . \end{array} . \end{aligned} \tag{5.89}$$

The components of T*<sup>R</sup>* are therefore given by

$$\begin{split} T\_{R\_{11}} &= \quad \left(1 + (c\_R - c\_0)\Omega\right) \left(\frac{\partial u\_R}{\partial R} + 1\right)^{-1} \left[2G\ln\left(\frac{\frac{\partial u\_R}{\partial R} + 1}{\left(1 + (c\_R - c\_0)\Omega\right)^{\frac{1}{3}}}\right) \right. \\ &\left. + \left(K - \frac{2}{3}G\right) \left(\ln\left(\frac{\frac{\partial u\_R}{\partial R} + 1}{\left(1 + (c\_R - c\_0)\Omega\right)^{\frac{1}{3}}}\right) \right. \\ &\left. + 2\ln\left(\frac{\frac{u\_R}{R} + 1}{\left(1 + (c\_R - c\_0)\Omega\right)^{\frac{1}{3}}}\right) \right) \right], \end{split} \tag{5.90}$$

$$\begin{split} T\_{R\_{22}} &= \quad \left(1 + (c\_R - c\_0)\Omega\right) \left(\frac{u\_R}{R} + 1\right)^{-1} \left[2G\ln\left(\frac{\frac{u\_R}{R} + 1}{\left(1 + \left(c\_R - c\_0\right)\Omega\right)^{\frac{1}{3}}}\right) \right. \\ &\left. + \left(K - \frac{2}{3}G\right) \left(\ln\left(\frac{\frac{\partial u\_R}{\partial R} + 1}{\left(1 + \left(c\_R - c\_0\right)\Omega\right)^{\frac{1}{3}}}\right) \right. \\ &\left. + 2\ln\left(\frac{\frac{u\_R}{R} + 1}{\left(1 + \left(c\_R - c\_0\right)\Omega\right)^{\frac{1}{3}}}\right) \right) \right]. \end{split} \tag{5.91}$$

Now we derive the coupled diffusion equation of the spherically symmetric boundary value problem. First, the coupled chemical potential µ *cp R* (*cR*,E *e log*), as shown in Equation (5.67), is expressed as

$$\begin{split} \mu\_{\boldsymbol{\kappa}}^{\rm cp}(\boldsymbol{c}\_{R},\mathbf{E}\_{\boldsymbol{\kappa}\_{\rm log}}^{\epsilon}) &= \quad G\boldsymbol{\Omega} \left[ \left( \ln \lambda\_{1}^{\epsilon} \right)^{2} + 2 \left( \ln \lambda\_{2}^{\epsilon} \right)^{2} + \frac{2\nu}{1 - 2\nu} \left( \ln \lambda\_{1}^{\epsilon} + 2 \ln \lambda\_{2}^{\epsilon} \right)^{2} \right] \\ &\quad + G \left( 1 + \left( \boldsymbol{c}\_{R} - \boldsymbol{c}\_{0} \right) \boldsymbol{\Omega} \right) \left[ \frac{-2\boldsymbol{\Omega} \left( 1 + \frac{\partial \boldsymbol{u}\_{R}}{\partial R} \right) \ln \lambda\_{1}^{\epsilon}}{3 \left( 1 + \left( \boldsymbol{c}\_{R} - \boldsymbol{c}\_{0} \right) \boldsymbol{\Omega} \right)^{\frac{2}{3}} \lambda\_{1}^{\epsilon}} \\ &\quad - \frac{4\boldsymbol{\Omega} \left( 1 + \frac{\boldsymbol{u}\_{R}}{R} \right) \ln \lambda\_{2}^{\epsilon}}{3 \left( 1 + \left( \boldsymbol{c}\_{R} - \boldsymbol{c}\_{0} \right) \boldsymbol{\Omega} \right)^{\frac{4}{3}} \lambda\_{2}^{\epsilon}} + \frac{2\nu}{1 - 2\nu} \left( \ln \lambda\_{1}^{\epsilon} + 2 \ln \lambda\_{2}^{\epsilon} \right) \\ &\quad - \frac{\left( \boldsymbol{-$$

Substituting Equations (5.9), (5.23) and (5.92) into Equation (5.68), we can derive the dimensionless coupled diffusion equation

$$\begin{split} 0 &= \quad \bar{R}^2 \frac{\partial \bar{c}\_R}{\partial \bar{t}} + \frac{\partial}{\partial \bar{R}} \left[ \frac{\partial}{\partial \bar{R}} \left( -\bar{R}^2 M(c\_R) \left( \mu\_\kappa^{m\nu p} + \mu\_\kappa^{gd} + \mu\_\kappa^{cp} (c\_R, \mathbf{E}\_{\log}^{\epsilon}) \right) \right) \right] . \end{split} \tag{5.93}$$

For the detailed derivation of the final dimensionless coupled diffusion equation, please see Appendix A.3.2.

# **6 Phase-field modeling of sodium-ion battery particles**

In chapter 4, we introduced the phase field theory related to species insertion into or extraction from a cathodic particle where phase segregation possibly occurs during discharging or charging. In order to take the mechanical deformation into account, different elastic strain energy functionals were motivated from different mechanical theories in chapter 5. In this chapter, a phase-field model for NaFPO is studied for the first time (see also Zhang and Kamlah [89]). We implemented the model in COMSOL Multiphysics<sup>r</sup> for a spherically symmetric problem of sodium insertion into or extraction from a cathodic particle made NaFPO. We first describe the determination of the material parameters for NaFPO. The model captures the important feature of phase segregation into a sodium-poor phase *FePO*<sup>4</sup> and a sodium-rich phase *Na*2/3*FePO*4. For this material, there is a visible difference for the concentration and stress between SDT and the finite deformation theories. Furthermore, we compare the two cathode materials NaFPO and LiFPO of lithium-ion batteries to each other in terms of phase changes and stresses, and show that, for the same stage of average concentration, although the miscibility gap of NaFPO is smaller than that of LiFPO, the stresses in the cathode material NaFPO are higher in the phase segregated state. As a result, the suppression of phase segregation by the elastic strain energy is more easily achieved in NaFPO compared to LiFPO.

## **6.1 Material parameters of Na***x***FePO**<sup>4</sup>

We consider a typical active nanoparticle of size *R*<sup>0</sup> = 150 *nm* for the particle radius. According to Nakayama et al. [120], the diffusion coefficient of Naions in NaFPO is an order of magnitude lower than that of Li-ions in LiFPO, which also can be confirmed by comparison to Zhu et al.[36]. The diffusion coefficient of LiFPO is 1×10−<sup>14</sup> *m* <sup>2</sup>/*s* [119, 121], therefore, we use the value *D*<sup>0</sup> = 1×10−<sup>15</sup> *m* <sup>2</sup>/*s* as the diffusion coefficient of NaFPO. In the following, we will specify the values of the other parameters for NaFPO including *cmax*, α1, α2, ¯λ, Ω, and *E*0.

#### **6.1.1 Determination of c***max*

The maximum concentration *cmax* with units of *mol*/*m* <sup>3</sup> means the maximum content of a certain species the host material can accept, which can be expressed by [116]

$$c\_{\text{max}} = \frac{1}{V\_0 N\_A},\tag{6.1}$$

where *V*<sup>0</sup> is the volume occupied by one species atom. This can be obtained by

$$V\_0 = \frac{V\_{cell}}{n},\tag{6.2}$$

with *n* being the number of atoms of the species that the host material can accept in a unit cell of volume *Vcell*.

According to Ref. [122], the structures of olivine LiFPO and olivine NaFPO are the same. Based on the crystal structure of LiFPO as shown in Ref. [123] and the lattice parameters in Table 6.1, we can calculate the maximum lithium concentration of LiFPO by Equation (6.1), which gives 2.29 × 10<sup>4</sup> *mol*/*m* 3 , matching the value used by Zeng and Bazant [119], and Delacourt and Safari [124]. Therefore, in the same way, the maximum sodium concentration of


NaFPO can be obtained as *cmax* = 2.1 × 10<sup>4</sup> *mol*/*m* 3 , which is smaller than that of LiFPO due to the larger unit cell volume *Vcell* for NaFPO.

Table 6.1: Lattice parameters of LiFPO and NaFPO.

### **6.1.2 Determination of** α<sup>1</sup> **and** α<sup>2</sup>

According to Lu et al. [37], at room temperature, the phase diagram of olivine NaFPO consists of two regions. For 0 < *x* < 2/3, phase segregation into a sodium-poor phase *FePO*<sup>4</sup> and a sodium-rich phase *Na*2/3*FePO*<sup>4</sup> is found to be favorable, which means that this is a two-phase region. (Following common praxis, we replace in chemical formulas dimensionless concentration ¯*c* by *x*.) Based on the unit cell parameters and volume as functions of *x* in NaFPO as shown in Ref. [37], two different values for the unit cell parameters and volume, respectively, can be identified when *x* increases to around 0.08, which is the manifestation of the occurrence of the phase segregation. Therefore, we can obtain the normalized sodium concentration for the initiation of phase segregation, which is around 0.08. In contrast to LiFPO, where transformation from *FePO*<sup>4</sup> into *LiFePO*<sup>4</sup> occurs directly, the system of NaFPO goes through an intermediate state at *Na*2/3*FePO*4. For 2/3 < *x* < 1 there is the solid-solution phase *NaxFePO*<sup>4</sup> which is a single-phase region. With regard to the homogeneous free energy density in the two-phase region of NaFPO, we assume

$$\begin{split} \left| \Psi^{m \ast p} (\vec{c}) \right| &= \left| k\_B T\_{ref} N\_{A} c\_{max} \left( \alpha\_1 \vec{c} + \frac{\alpha\_2}{2} \vec{c}^2 \right. \\ &\left. + \frac{T}{T\_{ref}} \left( \vec{c} \ln \vec{c} + \left( \frac{2}{3} - \vec{c} \right) \ln \left( \frac{2}{3} - \vec{c} \right) \right) \right) . \end{split} \tag{6.3}$$

This is the classical homogeneous free energy density function for a two phase material [27, 85], which we have formulated such that it is limited to the range 0 < *x* < 2/3. Consequently, we fit the homogeneous free energy density (6.3) to the above experimental data [37] with respect to the parameters α<sup>1</sup> and α2, as shown in Fig. 6.1. The best fit was achieved with α<sup>1</sup> = 5 and α<sup>2</sup> = −15 at room temperature.

Figure 6.1: Fit of the normalized homogeneous free energy density in the two-phase region to the experimental data from Ref. [37] [89].

Fig. 6.1 shows the fitting result with the normalized homogeneous free energy density in the two-phase region ψ¯ *mwp* plotted versus normalized concentration *c*¯. It can be seen that ψ¯ *mwp* exhibits a doublewell structure, such that two different relative minima A and B occur, characterizing two phases of different solubility limits *c*<sup>α</sup> and *c*<sup>β</sup> respectively, of the sodium concentration. The two minima A and B correspond to the two phases *FePO*<sup>4</sup> and *Na*2/3*FePO*4, respectively. The two zones AC and BD are the "nucleation zones", and phase segregation occurs upon sufficient disturbance of the system's equilibrium. In the inner zone of concavity between the two points of inflection C and D, which is denoted the "spinodal decomposition zone", homogeneous sodium concentration states are unstable and phase segregation is initiated in any case. The Maxwell construction represents the volume fractions of a phase segregated system. The inflection points correspond to the second order derivative of ψ¯ *mwp* with respect to ¯*c*, as shown in Fig. 6.2, and the inflection point C matches the experimental data from Ref. [37].

Figure 6.2: Second order derivative of ψ¯ *mwp* with respect to ¯*c* [89].

### **6.1.3 Determination of** λ

The so-called reference value λ<sup>0</sup> for λ can be estimated as [116]

$$
\lambda\_0 = \frac{-\alpha\_2 \left(\frac{1}{N\_A c\_{\text{max}}}\right)^{\frac{2}{3}}}{3},
\tag{6.4}
$$

and the relation between λ<sup>0</sup> and λ is

$$
\lambda = \alpha \lambda\_0,\tag{6.5}
$$

where α is constant.

Using Equation (6.4), we obtain λ<sup>0</sup> = 9.1×10−<sup>19</sup> *m* 2 for NaFPO. In absence of any further information, we pick α = 20 as the value for NaFPO which is the one digit number in the order of magnitude of the values for LiFPO and LMO, see Table 6.2. Therefore, we get λ = 1.8 × 10−<sup>17</sup> *m* 2 for NaFPO. The values of λ for three cathode materials are shown in Table 6.2, as well.


Table 6.2: Gradient energy coefficients of three cathode materials.

#### **6.1.4 Determination of** Ω

The partial molar volume Ω plays a role analogous to a thermal expansion coefficient, meaning we can calculate Ω by the relation ε *<sup>s</sup>* = 1/3Ω(*c* − *c*0)I. Thus, based on Table 6.3, the partial molar volume of LiFPO is obtained as

$$
\Omega = \frac{0.022 \times 3}{(1 - 0) \, c\_{\text{max}}} = 2.9 \times 10^{-6} \, m^3/mol,\tag{6.6}
$$

which matches the value used by Song et al. [117]. Therefore, in the same way, the partial molar volume of NaFPO can be calculated as

$$
\Omega = \frac{0.041 \times 3}{\left(\frac{2}{3} - 0\right) c\_{\max}} = 8.8 \times 10^{-6} \text{ m}^3/mol. \tag{6.7}
$$

This value of Ω is larger than that of LiFPO, which is consistent with the fact that sodium has a larger cation radius than lithium as shown in Table 1.1. It should be noticed that we do not consider the possibility of a concentration dependence of Ω.


Table 6.3: Volume change of LiFPO and NaFPO.

#### **6.1.5 Determination of E**<sup>0</sup>

Young's modulus of the host material *FePO*<sup>4</sup> is 125 *GPa*, Young's modulus of *LiFePO*<sup>4</sup> is 123.9 *GPa*, and Young's modulus of the cathode material LiFPO is 124.5 *GPa*, which is the average of the values for *FePO*<sup>4</sup> and *LiFePO*<sup>4</sup> [126]. Based on the above data, due to the larger ionic radius of the Na-ion compared to that of a Li-ion, we regard the two digit estimate 120 *GPa* as Young's modulus of NaFPO. The material parameters for the two cathode materials LiFPO and NaFPO are summarized in Table 6.4.


Table 6.4: The material parameters for the cathode materials LiFPO and NaFPO.

## **6.2 Numerical implementation**

The resulting set of equations has been implemented in COMSOL Multiphysics<sup>r</sup> for solution by the finite element method. For the numerical methods, since the Cahn-Hilliard equation is a fourth order partial differential equation for the species concentration field, the standard finite element method with C<sup>0</sup> -continuous Lagrange basis functions is not sufficient for discretization. There are several methods to overcome these numerical difficulties. For a straightforward treatment of the high-order operator, Zhao et al. [97] employed the isogeometric analysis, which allows for the employment of a wide range of smooth, higher-order basis functions, providing global C<sup>1</sup> -continuity. An alternative approach is the use of split-methods or mixed-methods to reduce the fourth order equation into two second-order equations. For example, Miehe et al. [127] considered the chemical potential as an additional degree of freedom to derive the mixed form for C<sup>0</sup> -continuous basis functions. Huttin and Kamlah [27] introduced the second derivative of the species concentration as unknown function to split the Cahn-Hilliard equation into two second-order equations, avoiding the problem with the computation of the logarithmic entropic terms in the homogeneous free energy density as shown in Equation (4.4). In addition, nonlocal species concentration [128] or micromorphic species concentration [94] is introduced as an extra degree of freedom to obtain an efficient implementation with C<sup>0</sup> -continuous Lagrange basis functions. In our implementations, we use the second derivatives of the species concentration as an additional degree of freedom.

## **6.3 Cahn-Hilliard model without mechanics**

We consider the quasistatic insertion and extraction of sodium for a particle of NaFPO at *C* = 0.001, see the boundary condition (5.30). With this C-rate, we study the behavior for dynamic, i.e. continuous insertion very close to a sequence of equilibrium states. In this way, the system is allowed to move along a path of relaxed quasiequilibrium states. In our simulations, we exclusively focus on the two-phase region of NaFPO (0 < *x* < 2/3). The average concentration *cavg*, also called "state of charge" (SOC) is defined as *cavg* = R <sup>B</sup> *cdV* ¯ /*V*.

Figure 6.3: Normalized average system free energy Ψ¯ *avg* and normalized homogeneous free energy density ψ¯ *mwp* as function of *cavg* and ¯*c*, respectively [89].

In Fig. 6.3, the markers represent values of the normalized average system free energy

$$
\Psi\_{avg} = \frac{1}{V} \int\_{\beta \mathcal{\mathbb{B}}} (\Psi^{m \nu p} + \Psi^{gd}(\text{grad } \bar{c}))dV \tag{6.8}
$$

during insertion and extraction, respectively, as function of *cavg*. For demonstration purposes, the plot of the normalized homogeneous free energy density vs. dimensionless concentration is entered, as well. Markers on the normalized homogeneous free energy density curve correspond to homogeneous states whereas markers nearby the path of the Maxwell construction correspond to phase segregated states, as they are illustrated in Figs. 6.4 and 6.5. Fig. 6.4

Figure 6.4: Normalized sodium concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> for different time instants during sodium insertion [89].

shows the normalized sodium concentration along the radial coordinate at different time instants during sodium insertion. We find that at the beginning of sodium insertion the concentration is homogeneous and the corresponding concentration is 0. Once *cavg* gets close to 0.08, which corresponds to the inflection point C in Fig. 6.1, phase segregation is initiated. Two different phases can be recognized in the particle, namely the sodium-poor phase *FePO*<sup>4</sup> in the center corresponding to the first minimum A and the sodium-rich phase *Na*2/3*FePO*<sup>4</sup> at the outside corresponding to the second minimum B. A smooth but very narrow interface with concentration changing rapidly but continuously separates them. When *cavg* grows up to 2/3, the intermediate phase *Na*2/3*FePO*<sup>4</sup> occupies all of the particle.

Figure 6.5: Normalized sodium concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> for different time instants during sodium extraction [89].

Next, we consider the sodium extraction case shown in Fig. 6.5. At the beginning the system is in a homogeneous state and the corresponding concentration is 2/3. Once *cavg* is reduced to around 0.6, which corresponds to the inflection point D in Fig. 6.1, phase segregation is initiated during sodium extraction. In contrast to the insertion case, the sodium-poor phase *FePO*<sup>4</sup> is at the outside while the sodium-rich phase *Na*2/3*FePO*<sup>4</sup> is in the center. At the end of sodium extraction, the sodium-rich phase *Na*2/3*FePO*<sup>4</sup> vanishes. It should be noticed that although to the particle surface a mass flux is applied, the concentration on the particle surface nearly stays at a constant value. This is due to the fact that because of the extremely low *C*-rate of *C* = 0.001 for quasistatic insertion and extraction of sodium, the system moves along a path of relaxed quasiequilibrium states.

### **6.4 Cahn-Hilliard model with mechanics**

Concerning mechanics, we first study the coupling between the Cahn-Hilliard theory and SDT, and then compare the different small and finite deformation theories. In order to study the coupling effect through the coupling energy on sodium diffusion and stress in the particle, we introduce a normalized Young's modulus as

$$
\bar{E} = \frac{E}{E\_0},
\tag{6.9}
$$

where *E*<sup>0</sup> is the value of Young's modulus shown in Table 6.4. In the figures, *c*, *c<sup>R</sup>* denote concentration of small and finite deformation elasticity, respectively, and *T<sup>H</sup>* = 1/3*Tii* is the hydrostatic stress in terms of the Cauchy stress.

### **6.4.1 Cahn-Hilliard model with small deformation theory**

Figs. 6.6 and 6.7 show normalized sodium concentration and hydrostatic stress along the radial coordinate at a stage of average concentration of *cavg* = 0.5 during sodium insertion based on SDT for different *E*¯. We find that the difference in concentration between the two phases is reduced as *E*¯ increases, which means that the solid solution limits of *FePO*<sup>4</sup> and *NaxFePO*<sup>4</sup> are gradually extended into the range of phase segregated states. In other words the coupling effect reduces the miscibility gap. Correspondingly, the hydrostatic stress magnitudes increase as *E*¯ increases. The hydrostatic stress is constant and tensile in the low concentration phase, and it drops rapidly across the interface to become compressive and constant in the high concentration phase. This change of sign is a result of mechanical equilibrium between the inner core and the outer shell of the particle. Note that on the other side for *E*¯ = 1,

the system is homogeneous and stress-free. Even though the average concentration value lies in the spinodal decomposition zone of ψ¯ *mwp*, phase segregation does not occur.

Figure 6.6: Normalized sodium concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 based on SDT for different *E*¯ during sodium insertion [89].

Figure 6.7: Hydrostatic stress *T<sup>H</sup>* versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 based on SDT for different *E*¯ during sodium insertion [89].

As illustrated in Fig. 6.8, for a high value of normalized Young's modulus, the cost of the system coupling energy Ψ¯ *cp* at a phase segregated state would be too high so that as a result of the competition between the different contributions to the system free energy the system stays in the homogeneous state in order to minimize the total system free energy. Therefore, the presence of a coupling energy in the total system free energy can lead to suppression of phase segregation. For the detailed discussion of the influence of mechanics on the miscibility gap, see chapter 7.

Figure 6.8: Normalized system coupling energy Ψ¯ *cp* versus average concentration *cavg* based on SDT for different *E*¯ during sodium insertion. Here Ψ¯ *cp* = R <sup>B</sup> ψ¯ *cpdV* [89].

#### **6.4.2 Comparison of different mechanics theories**

Figs. 6.9 and 6.10 show the comparison of normalized sodium concentration and hydrostatic stress at a stage of average concentration of *cavg* = 0.5 during sodium insertion for different mechanics theories with a value of *E*¯ = 0.3: SDT, elasticity based on elastic Green strain, and elasticity based on logarithmic elastic strain. According to the experimental work of Xiang et al. [39], driven by the coupling energy, the sodium-rich phase decreases its sodium content to about *Na*0.6*FePO*4, while the sodium-poor phase increases its sodium content to about *Na*0.08*FePO*4. These experimental results [39] are shown as dashed horizontal lines in Fig. 6.9. It can be recognized that there is a visible difference for the sodium concentration plots of SDT and the finite deformation elasticity formulations, with the difference between the two phases being less pronounced for SDT, while the sodium concentration plots of two finite deformation theories are almost the same. Compared to the experimental result [39], we can find that the sodium concentration in the high concentration phase obtained from two finite deformation theories matches the experimental value for *Na*0.6*FePO*<sup>4</sup> but there is a clear deviation from the result by SDT. On the other hand, the sodium concentration in the low concentration phase obtained from all three mechanics theories is always slightly below the experimental value for *Na*0.08*FePO*4. In Fig. 6.10, the hydrostatic stress plots are clearly different for SDT on the one side and the two finite deformation theories on the other. There are higher stresses in the two phases for SDT. This is attributed to the fact that SDT does not account for the volume swelling of the particle during sodium insertion. Comparing the two finite deformation theories to each other, it is found that there is a negligible difference in stress between them.

Figure 6.9: Normalized sodium concentration versus normalized radial coordinate at *cavg* = 0.5 with *E*¯ = 0.3 during sodium insertion for different mechanics theories. The normalized radial coordinates *r*/*R*<sup>0</sup> and *R*/*R*<sup>0</sup> represent SDT and finite deformation theory, respectively. The dashed horizontal lines represent the experimental results from Ref. [39] [89].

As mentioned before, the contribution of the coupling energy to the total system free energy can lead to suppression of phase segregation. Here, we discuss the critical value of normalized Young's modulus *E*¯ *<sup>c</sup>* of NaFPO above which phase segregation can not arise in the particle. Fig. 6.11 shows the normalized sodium concentration along the radial coordinate for various values of *E*¯ during sodium insertion based on different mechanics theories. For SDT, it can be seen that phase segregation can arise at a stage of average concentration of *cavg* = 0.33 when *E*¯ = 0.385 but sodium concentration is homogeneous over the particle when *E*¯ = 0.386, so *E*¯ *<sup>c</sup>* is 0.385 for SDT. In the same way, we find *E*¯ *<sup>c</sup>* = 0.409 for the two finite deformation elasticity formulations, which is a little larger than that for SDT.

Figure 6.10: Hydrostatic stress versus normalized radial coordinate at *cavg* = 0.5 with *E*¯ = 0.3 during sodium insertion for different mechanics theories [89].

## **6.5 Comparison with the cathode material Li***x***FePO**<sup>4</sup>

We now compare the two cathode materials NaFPO and LiFPO in terms of phase changes and hydrostatic stress. Figs. 6.12 and 6.13 show normalized concentration and hydrostatic stress along the radial coordinate for various particle radii *R*<sup>0</sup> for the two cathode materials during insertion based on different finite deformation theories. Here, we consider two stages of average concentration, one is *cavg* = (*c*<sup>α</sup> + *c*<sup>β</sup> )/2 which is the center of the combined nucleation and spinodal zone, the other is *cavg* = 0.5. It should be noticed that

Figure 6.11: Normalized sodium concentration versus normalized radial coordinate for various *E*¯ during sodium insertion based on different mechanics theories [89].

(*c*<sup>α</sup> + *c*<sup>β</sup> )/2 of LiFPO is just equal to 0.5, and which of NaFPO is 0.333. In Fig. 6.12, we recognize that there is an obviously reduced miscibility gap for a NaFPO particle compared to a LiFPO particle. This is consistent with the experimental fact that NaFPO has a two-phase region in the range 0 < *x* < 2/3 while phase segregation for LiFPO occurs in the range 0 < *x* < 1 [37]. For *cavg* = (*c*<sup>α</sup> +*c*<sup>β</sup> )/2, the two cathode materials show almost the same position of the center of the interface region, which is located close to the particle surface due to the spherical symmetry. Additionally, for *cavg* = 0.5, the interface location in a NaFPO particle is more close to the particle center, and, furthermore, the concentration in the two phases is the same for different particle radii. As shown in Fig. 6.13, for *cavg* = (*c*<sup>α</sup> + *c*<sup>β</sup> )/2, both the tensile stresses in the inner core and the compressive stress magnitudes in the outer shell of a NaFPO particle are larger than those in a LiFPO particle. When *cavg* increases from 0.333 to 0.5, the tensile stresses in the inner core of a NaFPO particle become larger, while the compressive stress magnitudes in the outer shell of a NaFPO particle decrease, which is even smaller that those in a LiFPO particle. Indeed, the particle shell made of high-concentration phase becomes thicker when *cavg* shifts from 0.333 to 0.5 (please see Fig. 6.12), so the particle core made of low-concentration phase is in mechanical equilibrium with a thicker shell made of high-concentration phase. As a result, both the tensile stresses in the inner core and the compressive stresses in the outer shell of a NaFPO particle increase. Therefore, for the same stage of average concentration, although the miscibility gap of NaFPO is much smaller than that of LiFPO, the stresses in a NaFPO particle are larger compared to those in a LiFPO particle. This is due to a larger expansion during phase changes for NaFPO as shown in Fig. 6.14 describing the plots of volume ratio *J* at the particle surface. Here, it should be mentioned again that our simulations only consider the two-phase region of NaFPO (0 < *x* < 2/3). We can find that the volume expansion from *FePO*<sup>4</sup> to *Na*2/3*FePO*<sup>4</sup> is quite large, namely about 12.8%, which is nearly 2 times that for LiFPO changing from *FePO*<sup>4</sup> to *LiFePO*4. This is consistent with the reports [37–39]. Correspondingly, as shown in Fig. 6.15, the maximum hydrostatic stress magnitude in a NaFPO particle during phase changes is always larger than that in a LiFPO particle. The larger stresses in a NaFPO particle may explain the existence of a wide range of solid solution *NaxFePO*<sup>4</sup> (2/3 < *x* < 1) to avoid even higher stresses in the NaFPO particle. Casas-Cabanas et al. [38] conclude that the larger stresses in a NaFPO particle can be the explanation for the existence of an intermediate phase. Indeed, the formation of an intermediate phase acts as a buffer between *FePO*<sup>4</sup> and *NaFePO*<sup>4</sup> providing elasticity to the structure. On the other hand, the critical value of normalized Young's modulus *E*¯ *<sup>c</sup>* of LiFPO for surpression of phase segregation is equal to 1 (please see Fig. 7.5 in chapter 7), which is larger than that of NaFPO, as shown before. Therefore, it is easier for NaFPO to reach a massive coupling energy to suppress phase segregation compared to LiFPO due to the larger volume change of a NaFPO particle.

As can be seen in Fig. 6.15, the LiFPO particle surface expands more seriously around *cavg* = 0.15. This is due to phase segregation being already initiated for LiFPO but not yet for NaFPO which is still in the homogeneous state. As a result, the maximum hydrostatic stress magnitude in a LiFPO particle is far greater than that in a NaFPO particle at this insertion state. At *cavg* = 2/3, the maximum hydrostatic stress magnitude in a NaFPO particle is close to zero while the LiFPO particle still displays high stresses. This is attributed to the intermediate phase *Na*2/3*FePO*<sup>4</sup> occupying then the whole particle while LiFPO is still in a phase segregated state.

Figure 6.12: Normalized concentration ¯*c<sup>R</sup>* versus normalized radial coordinate *R*/*R*<sup>0</sup> with *E*¯ = 0.3 for various *R*<sup>0</sup> during insertion for two cathode materials based on different finite deformation elasticity formulations. Two stages of average concentration are considered: *cavg* = (*c*<sup>α</sup> +*c*<sup>β</sup> )/2 and *cavg* = 0.5, where (*c*<sup>α</sup> +*c*<sup>β</sup> )/2 of LiFPO is just equal to 0.5, and which of NaFPO is 0.333.

Figure 6.13: Hydrostatic stress *T<sup>H</sup>* versus normalized radial coordinate *R*/*R*<sup>0</sup> with *E*¯ = 0.3 for various *R*<sup>0</sup> during insertion for two cathode materials based on different finite deformation elasticity formulations. Two stages of average concentration are considered: *cavg* = (*c*<sup>α</sup> +*c*<sup>β</sup> )/2 and *cavg* = 0.5, where (*c*<sup>α</sup> +*c*<sup>β</sup> )/2 of LiFPO is just equal to 0.5, and which of NaFPO is 0.333.

Figure 6.14: Volume ratio *J* at the particle surface versus *cavg* with *E*¯ = 0.3 for various *R*<sup>0</sup> during insertion for two cathode materials based on different finite deformation elasticity formulations [89].

Figure 6.15: Maximum hydrostatic stress magnitude |*TH*,*max*| versus *cavg* with *E*¯ = 0.3 for various *R*<sup>0</sup> during insertion for two cathode materials based on different finite deformation elasticity formulations [89].

# **7 Particle size and average concentration dependent miscibility gap**

In the previous chapter 6, we just focus on the bulk particle size of the cathode material. As mentioned before, the thermodynamics of phase segregation in nanoscale particles is distinctly different from bulk materials, and the particle size miscibility gap plays a nontrivial role in the performance of nanoscale insertion materials. In contrast to our common thermodynamic knowledge, the miscibility gap in nanoscale insertion material also depends on the average concentration. In this chapter, the particle size and average concentration dependent miscibility gap of the nanoscale insertion materials LMO, LiFPO, and NaPPO are investigated during insertion, using a coupled phase-field model based on the Cahn-Hilliard theory and finite deformation elasticity (see also Zhang and Kamlah [96]). For the mechanical part, two finite deformation elasticity formulations, based on elastic Green strain and logarithmic elastic strain, respectively, have been compared in the previous chapter 6. We found that, although the results from the two formulations of Hooke's law with the geometric nonlinearity of finite strain are almost the same, according to our experience, elasticity based on logarithmic elastic strain is numerically more efficient than elasticity based on elastic Green strain. The stress strain relation for LiFPO is also described by elasticity based on logarithmic elastic strain in the work by Di Leo et al. [94]. Therefore, a logarithmic elastic strain law is employed in this chapter. We implemented the model in COMSOL Multiphysics<sup>r</sup> for a spherically symmetric problem of species insertion into cathodic particles which are respectively made of LMO, LiFPO, and NaPPO. We first compare three bulk cathode materials in terms of phase changes and stresses, and then investigate their particle size and average concentration dependent miscibility gap.

### **7.1 Material parameters**

In order to study the particle size and average concentration dependent miscibility gap in different nanoinsertion cathode materials, the three cathode materials LMO, LiFPO, and NaFPO are chosen in the simulations. The material parameters for these three cathode materials are summarized in Table 7.1, where the material parameters of NaFPO have been determined in chapter 6. We vary the nanoparticle radius to study the influence of particle size on the miscibility gap. Depending on the temperature, LMO possibly exhibits phase segregation in the range of values 0 < *x* < 1 where the crystalline host structure remains cubic spinel [132–134]. Some theoretical calculations of LMO have been carried out based on an assumption that Li-ordering phase exists at *x* = 0.5 around room temperature [23, 135, 136], but this does not agree with the experimental literature that suggest no Li-ordering phase [137]. At *x* = 1, LMO undergoes a phase transition where the crystalline material becomes tetragonal due to the Jahn-Teller distorsion regarding the manganese ions [138]. For LMO, we focus on the phase segregation that happens at room temperature in the range of values 0 < *x* < 1, neglecting the effect of lithium ordering. The cubic to tetragonal transition at *x* = 1 is also not taken into account since the state of charge is rarely driven up to *x* = 1 unless delivering a large current density [139].

At room temperature, the miscibility gap of bulk LiFPO encompasses nearly the entire lithium composition range of the material [50] compared to a relatively small miscibility gap of bulk LMO bounded by solid solutions states between *Li*0.12*Mn*2*O*<sup>4</sup> and *Li*0.88*Mn*2*O*<sup>4</sup> [116]. However, bulk NaFPO shows


Table 7.1: The material parameters for three cathode materials.

the smallest miscibility gap bounded by solid solutions states between *FePO*<sup>4</sup> and *Na*2/3*FePO*<sup>4</sup> [37] among these three cathode materials.

The three cathode materials show different volume changes during species insertion. For LMO, the volume expansion of *Mn*2*O*<sup>4</sup> upon lithiation to *LiMn*2*O*<sup>4</sup> reaches about 7.3% [140]. According to Lu et al. [37], the volume expansion of *FePO*<sup>4</sup> upon sodiation to *NaFePO*<sup>4</sup> is quite large, namely about 17%, and even the volume expansion from *FePO*<sup>4</sup> to *Na*2/3*FePO*<sup>4</sup> is still quite large (about 12.8%), which is nearly 2 times that for LiFPO (about 6.8%) changing from *FePO*<sup>4</sup> to *LiFePO*4. The previous chapter 6 reveals that the difference for the concentration and stress plots between SDT and finite deformation elasticity can not be neglected for NaFPO. Therefore, finite deformation elasticity is selected in order to represent all the three cathode materials by the same deformation theory.

## **7.2 Comparison of the three bulk cathode materials**

We will consider quasistatic insertion of species into nanoparticles for three cathode materials at *C* = 0.001. Again, Ψ¯ *avg* = R B<sup>R</sup> ψ¯*RdVR*/*V<sup>R</sup>* is the normalized average system free energy, see also Equation (5.48). We compare three bulk cathode materials in terms of phase changes and stresses, and a typical active particle radius of *R*<sup>0</sup> = 200 *nm* is considered.

Fig. 7.1 shows the values of the normalized average system free energy for the three cathode materials during insertion by the markers, including the pure diffusion case without mechanical deformation and the mechanically coupled diffusion case, as function of *cavg*. For demonstration purposes, the plots of the normalized homogeneous free energy density versus the normalized concentration are also entered as line plots. The species-poor phase and the species-rich phase are characterized by the solubility limits *c*<sup>α</sup> and *c*<sup>β</sup> , respectively, of concentration. Here, the solubility limits are the lower and upper bounds of the concentration range of phase segregated states. For example, in Fig. 7.1, the solubility limits can be identified where the normalized average system free energy states do not coincide with the plot of the homogenous multiwell potential. The difference in the value of ψ¯ *mwp* at a stage of average concentration of *cavg* = (*c*<sup>α</sup> + *c*<sup>β</sup> )/2, and the value of ψ¯ *mwp* at either of the solubility limits *c*<sup>α</sup> or *c*<sup>β</sup> represents the maximum energy barrier ∆ψ¯ *mwp max* . As mentioned in section 6.1.2, the two ranges between the respective point of tangency defined by the Maxwell construction and the corresponding inflection point are the "nucleation zones", and the inner zone of concavity between the two points of inflection is called the "spinodal decomposition zone". Again, markers on the normalized homogeneous free energy density curve correspond to homogeneous states whereas markers nearby the path of the Maxwell construction correspond to phase segregated states. In particular, at a stage of average concentration of *cavg* = 0.5, all three cathode materials show phase segregated state, as illustrated in Fig. 7.2 (a).

Figure 7.1: Normalized average system free energy Ψ¯ *avg* and normalized homogeneous free energy density ψ¯ *mwp* as function of *cavg* and ¯*c*, respectively, for the three cathode materials and a radius of *R*<sup>0</sup> = 200 *nm*. The dashed straight lines represent Maxwell construction, and *E*¯ = 0.3 is considered in the mechanically coupled diffusion case [96].

We first focus on the pure diffusion case. In Fig. 7.1, we find that phase segregation is initiated for NaFPO once *cavg* approaches the inflection point at *cavg* = 0.08, which is earlier compared to LiFPO (*cavg* = 0.13) and LMO (*cavg* = 0.26). Also, it can be seen that LiFPO has the largest "spinodal decomposition zone" among the three cathode materials. In Fig. 7.2 (a), we also find that the miscibility gap, when neglecting mechanics, is the smallest for NaFPO, which is about 2/3 of that of LiFPO. However, if the contribution from the coupling energy is considered (*E*¯ = 0.3), the miscibility gap of LMO is significantly reduced. We can see that in this case the related normalized average system free energy shown in Fig. 7.1 coincides with the respective normalized homogeneous free energy density even in most of the "spinodal decomposition zone". Hardly visible, phase segregation, i.e. a minor deviation from ψ¯ *mwp* takes place in the neighborhood of *cavg* ≈ 0.5, only. The resulting miscibility gap is even smaller than that of NaFPO, even if mechanics is accounted for the latter material. The physical explanation for this behavior of LMO is that the cost of the elastic strain energy at a phase segregated state would be too high so that as a result of the competition between the different contributions to the system free energy the system stays in the homogeneous state in order to minimize the total system free energy. The effect of mechanics, i.e. the system coupling free energy, is the weakest for LiFPO. Indeed, as shown in Fig. 7.1 (a), in contrast to the other two materials, when considering mechanics, the markers of LiFPO are still close to the path of the Maxwell construction, and the miscibility gap of LiFPO is just slightly reduced in Fig. 7.2 (a). As a result of the larger miscibility gap, which by the partial molar volume translates into a larger strain mismatch, the tensile stresses in the inner core of a LiFPO particle are larger compared to those in a LMO particle, and the compressive stress magnitudes in the outer shell of a LiFPO particle are also larger than those in a LMO particle, as shown in Fig. 7.2 (b). On the other hand, as mentioned in section 6.5, although the miscibility gap of NaFPO is smaller than that of LiFPO, both the tensile stresses in the inner core and the compressive stresses in the outer shell of a NaFPO particle are larger than those in a LiFPO particle at *cavg* = 0.5. This again is due to a larger expansion, i.e. partial molar volume, during phase segregation for NaFPO.

Figure 7.2: (a) Normalized concentration versus normalized radial coordinate at *cavg* = 0.5 during insertion for three cathode materials of *R*<sup>0</sup> = 200 *nm*. The solid lines represent the pure diffusion case and the dashed lines represent the mechanically coupled diffusion case of *E*¯ = 0.3. (b) Hydrostatic stress *T<sup>H</sup>* versus normalized radial coordinate *R*/*R*<sup>0</sup> at *cavg* = 0.5 during insertion for three cathode materials (*R*<sup>0</sup> = 200 *nm* and *E*¯ = 0.3) [96].

### **7.3 Particle size dependent miscibility gap**

In order to investigate the impact of the particle size on the miscibility gap for the three cathode materials, we first exclude the effect of the mechanics. The solubility limits of the material are determined by taking the minimum and maximum concentrations for which phase segregated states occur [50].

Fig. 7.3 shows the solubility limits as a function of particle radius at a stage of average concentration of *cavg* = (*c*<sup>α</sup> +*c*<sup>β</sup> )/2 for three cathode materials. The motivation for focusing on this specific value of average concentration is based on observations in our simulation results, that the miscibility gap shrinks symmetrically about the center of the combined nucleation and spinodal zone.

Consequently, we assume that phase segregation is suppressed completely,

when it is suppressed for *cavg* = (*c*<sup>α</sup> + *c*<sup>β</sup> )/2. Here, *c*<sup>α</sup> and *c*<sup>β</sup> are the size independent solubility limits found for sufficiently large particles.

We find that a reduction of the particle radius leads to a shrinking of the miscibility gap. In particular, for each material a specific critical particle radius is found below which phase segregation is inhibited. The critical radius of LMO (*R<sup>c</sup>* = 10.9 *nm*) is larger than those of 6 *nm* for LiFPO and 6.4 *nm* for NaFPO, respectively. The above predicted critical particle sizes match the fact that phase segregation would be suppressed when the particle size becomes comparable to the interface thickness. According to Cahn and Hilliard [29], the interface thickness *d* is proportional to the square root of the gradient coefficient λ and the inverse of the maximum energy barrier , i.e *d* ∝ (λ/∆ψ¯ *mwp max* ) 1 2 . Therefore, LMO has the largest critical particle size among the three cathode materials due to a very tiny ∆ψ¯ *mwp max* shown in Fig. 7.1. It should be mentioned that our value of *R<sup>c</sup>* = 6 *nm* for LiFPO is a little bit smaller than the estimate value *R<sup>c</sup>* = 7.5 *nm* from Meethong et al. [45], but it is close to the critical radius (around 5 *nm*) obtained from a 3D comprehensive phase-field model accounting for facet dependent surface wetting [51]. On the other hand, different materials show different particle size dependent miscibility gap behavior. For LMO a reduced miscibility gap is found below a particle radius of 55.5 *nm*. The particle size dependent miscibility gap becomes more and more obvious as the particle size decreases, and strongly varying solubility limits are found below a particle radius of 35 *nm*. Compared to LMO, the threshold values of the particle radius for the size independent miscibility gap are smaller for LiFPO and NaFPO, which are 20.5 *nm* and 17.5 *nm*, respectively. Therefore, LMO has the broadest particle size range for varying solubility limits.

Figure 7.3: Solubility limits as a function of particle radius at *cavg* = (*c*<sup>α</sup> +*c*<sup>β</sup> )/2 during insertion for three cathode materials [96].

With the help of Fig. 7.4, we take LiFPO as an example to explain the shrinking of the miscibility gap upon particle size reduction. As the particle size reduces, the interface region, the thickness of which is determined by the material constant λ, covers a larger fraction of the material. Furthermore, due to the spherical symmetry, the position of the center of the interface region for *cavg* = (*c*<sup>α</sup> +*c*<sup>β</sup> )/2 ≈ 0.5 is located close to the particle surface for all values of radius. Consequently, below a radius of 20.5 *nm* the concentration value at the surface progressively is reduced from an initial value of *c*<sup>β</sup> . As a result the miscibility gap shrinks at the same average concentration to accommodate in a particle of reduced size the concentration gradient, the slope of which is basically given as a material constant. Due to the normalization of the radial coordinate in the form *r*/*R*<sup>0</sup> in Fig. 7.4, this effect translates into an apparent decrease of the concentration gradient, resulting in the same decrease of the surface concentration. Once the particles reduces to be below the critical particle size, the cost of the gradient energy Ψ¯ *gd* at a phase segregated state would be too high so that the system would stay in a homogeneous state in order to minimize the total system free energy due to the competition between the different contributions to the system free energy, namely multiwell potential and gradient energy.

Figure 7.4: Normalized concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* ≈ 0.5 during insertion for different particle sizes of LiFPO [96].

We now investigate to which extent mechanics affects the particle size dependent miscibility gap for the example of based LiFPO with *E*¯ = 1 for *cavg* = (*c*<sup>α</sup> + *c*<sup>β</sup> )/2 ≈ 0.5. In Fig. 7.5, we can find that the threshold value of the particle radius for the size independent miscibility gap expands to 43 *nm* in the mechanically coupled diffusion case, and the critical particle radius below which phase segregation is inhibited increases to 9 *nm*. Since the strain mismatch between regions of different phase state induces the additional elastic coupling energy which scales with the width of the miscibility gap through the partial molar volume, minimization of the total energy leads to a reduction of the miscibility gap when mechanics is accounted for. The finding of a reduced miscibility gap in the presence of mechanics is also consistent with the distribution of hydrostatic stress in Fig. 7.6 (a). There, the gradient of hydrostatic stress induces a mechanical contribution to the driving force for diffusion which points to the particle center, i.e. from the high concentration phase towards the low concentration phase. As a result, the hydrostatic stress assists species transport towards the particle center and the concentration fields considering the mechanics are more diffuse than that of the pure diffusion case, as shown in Fig. 7.6 (b). Fig. 7.6 (b) confirms that for the same particle radius not only the miscibility gap gets smaller, but also the width of the interface is widened when accounting for mechanics. This means that a comparably larger interface size needs to be accommodated inside the particle. This leads to a larger threshold value of the particle radius for the size independent miscibility gap and it leads to a larger particle size below which no phase segregation is obtained at all. On the other hand, for the small particle size the difference in concentration distribution between the pure and mechanically coupled diffusion cases increases as the particle size decreases.

Figure 7.5: Solubility limits as a function of particle radius at *cavg* ≈ 0.5 during insertion for the pure and mechanically coupled diffusion cases of LiFPO. *E*¯ = 1 is considered in the mechanically coupled diffusion case [96].

Figure 7.6: (a) Hydrostatic stress *T<sup>H</sup>* versus normalized radial coordinate *R*/*R*<sup>0</sup> at *cavg* ≈ 0.5 during insertion for different particle sizes of LiFPO with *E*¯ = 1. (b) Normalized concentration versus normalized radial coordinate at *cavg* ≈ 0.5 during insertion for different particle sizes of LiFPO, where the solid lines represent the mechanically coupled diffusion case of *E*¯ = 1 and the dashed lines represent the pure diffusion case [96].

## **7.4 Evolution of the miscibility gap during insertion**

As discussed above the size dependence of the miscibility gap has been investigated so far at a stage of average concentration of *cavg* = (*c*<sup>α</sup> + *c*<sup>β</sup> )/2. This value is of particular importance for the critical particle radius below which phase segregation is inhibited. We now investigate how for a fixed particle size the miscibility gap evolves during the process of insertion. Again, we first exclude the effect of mechanics.

Figs. 7.7 - 7.9 show that for the smaller particle sizes the solubility limits significantly depend on the average concentration for all three cathode materials. As the particle size decreases, an extended solid solution range is observed even in the metastable nucleation and the unstable spinodal regions of average concentration. In general, the miscibility gap expands as the average concentration increases. For larger particles, the system eventually enters a constant miscibility gap state. Interestingly, for *R*<sup>0</sup> = 12.5 *nm* the miscibility gap of LMO in the whole range of phase segregated states is even smaller than that of NaFPO, while for larger particle sizes it is the other way around. This means that the suppression of phase segregation by reducing the particle size is more easily achieved in LMO. This is also consistent with the fact that, among the three materials, LMO has the largest critical radius below which phase segregation is inhibited. On the other hand, once the particle radius reaches about 125 *nm*, LMO reveals a miscibility gap which is constant in the whole range of phase segregated states, i.e. a miscibility gap indepentdent of average concentration. Compared to LMO, the threshold values of the particle radius for the miscibility gap to be independent of average concentration are smaller for LiFPO and NaFPO, namely 73 *nm* and 75 *nm*, respectively.

Figure 7.7: (a) Solubility limits as a function of average concentration *cavg* during insertion for different particle sizes of LMO. The red and black curves represent *c*<sup>α</sup> and *c*<sup>β</sup> , respectively. (b) Normalized concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> for different time instants during insertion for three different particle sizes of LMO [96].

Figure 7.8: (a) Solubility limits as a function of average concentration *cavg* during insertion for different particle sizes of LiFPO. The red and black curves represent *c*<sup>α</sup> and *c*<sup>β</sup> , respectively. (b) Normalized concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> for different time instants during insertion for three different particle sizes of LiFPO [96].

Figure 7.9: (a) Solubility limits as a function of average concentration *cavg* during insertion for different particle sizes of NaFPO. The red and black curves represent *c*<sup>α</sup> and *c*<sup>β</sup> , respectively. (b) Normalized concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> for different time instants during insertion for three different particle sizes of NaFPO [96].

The evolution of the miscibility gap with average concentration for the smaller particle sizes can be explained by the minimization of system free energy for the example of LiFPO, as shown in Fig. 7.10. For *R*<sup>0</sup> = 100 *nm*, the system reaches the maximum miscibility gap immediately once phase segregation is initiated (*cavg* = 0.13 for *R*<sup>0</sup> = 100 *nm*), such that the system exhibits a miscibility gap, which is independent of average concentration. For smaller particle sizes, for example, *R*<sup>0</sup> = 12.5 *nm*, the system shows a dynamically increasing miscibility gap behavior until *cavg* = 0.49. What is the reason behind that? First, the interface volume is the largest at the beginning of phase segregation due to the shrinking core dynamics, leading to a high cost of average gradient energy (Ψ¯ *gd avg* = R <sup>B</sup> 1/2 λ/*R* 2 0 | ~∇¯ *c*¯| 2 *dV*/*V*). Second, there exists a relatively large average gradient energy for smaller particle sizes, which is mainly controlled by the ratio λ/*R* 2 0 . Therefore, if the maximum miscibility gap immediately occurs once phase segregation is initiated for smaller particle sizes, the cost of the total system free energy would be even higher than the normalized chemical potential. In order to minimize the system free energy, the system shows a dynamically increasing miscibility gap behavior for the smaller particle sizes before the system reaches the maximum miscibility gap.

Figure 7.10: Normalized average system free energy Ψ¯ *avg* and normalized homogeneous free energy density ψ¯ *mwp* as function of *cavg* and ¯*c*, respectively, for different particle sizes of LiFPO during insertion.

Now we investigate how mechanics influences the average concentration dependent miscibility gap based on LiFPO. As expected, the miscibility gap of the mechanically coupled diffusion case shrinks due to the suppressing effect of the coupling energy, as shown in Fig. 7.11. When accounting for mechanics, the initiation of phase segregation is postponed and the miscibility gap only exists for a smaller range of the average concentration. For the smaller particle radius, for example, *R*<sup>0</sup> = 12.5 *nm*, the miscibility gap of the mechanically coupled diffusion case always expands in almost the whole range of phase segregated states, while the miscibility gap of the pure diffusion case becomes constant once *cavg* approaches 0.49. On the other hand, the threshold value for the miscibility gap to be completely independent of average concentration expands to 100 *nm* in the mechanically coupled diffusion case compared to 73 *nm*of the pure diffusion case shown in Fig. 7.8 a.

Figure 7.11: Solubility limits as a function of average concentration *cavg* during insertion for different particle sizes of LiFPO. The solid lines represent the mechanically coupled diffusion case of *E*¯ = 1 and the dashed lines represent the pure diffusion case [96].

# **8 Nonlocal species concentration model**

In chapter 4, a phase-field theory named NSC theory in terms of diffusion and phase changes has been derived. This theory incorporates two secondorder partial differential equations involving second-order spatial derivatives of species concentration and an additional variable called nonlocal species concentration. In this chapter, we implemented the NSC theory in COMSOL Multiphysics<sup>r</sup> for a spherically symmetric boundary value problem of lithium insertion into a LMO cathode material particle of a LIB, and present the results for spherical LMO electrode particles based on the NSC theory (see also Zhang and Kamlah [98]). We first introduce the spherically symmetric boundary value problem for this theory. Then, the material parameters controlling the interface are determined for LMO and the interface evolution is studied. Comparison to the Cahn-Hilliard theory shows that nonlocal species concentration theory is superior when simulating problems where the dimensions of the microstructure such as phase boundaries are of the same order of magnitude as the problem size. This is typically the case in nanosized particles of phase separating electrode materials. For example, the nonlocality of nonlocal species concentration theory turns out to make the interface of the local concentration field thinner than in Cahn-Hilliard theory. Finally, the influence of interface related parameters, including the penalty energy coefficient β and the characteristic interface length scale *l*, is investigated.

## **8.1 Spherically symmetric boundary value problem**

We consider a spherical particle of radius *R*<sup>0</sup> under spherically symmetric boundary conditions. Here, spherical coordinates are introduced and all fields are expressed as a function of time *t* and the radial coordinate 0 ≤ *r* ≤ *R*<sup>0</sup> in the form

$$c = c(r, t),\tag{8.1}$$

$$
\hat{c} = \hat{c}(r, t). \tag{8.2}
$$

Figure 8.1: Boundary and symmetry conditions for the spherically symmetric boundary value problem of the NSC theory.

For the spherically symmetric boundary value problem, based on Equation (4.54), the dimensionless diffusion equation for the NSC theory can be obtained as

$$\begin{split} 0 &= \quad \bar{r}^2 \frac{\partial \bar{c}}{\partial \bar{t}} + \frac{\partial}{\partial \bar{r}} \left[ -\bar{r}^2 (1 + \alpha\_2 \bar{c} (1 - \bar{c})) \frac{\partial \bar{c}}{\partial \bar{r}} - \bar{\beta} \bar{r}^2 \bar{c} (1 - \bar{c}) (\frac{\partial \bar{c}}{\partial \bar{r}} - \frac{\partial \bar{\hat{c}}}{\partial \bar{r}}) \right], \end{split} \tag{8.3}$$

where the nonlocal species concentration is calculated by solving the dimensionless Helmholtz equation

$$\frac{\lambda}{R\_0^2} \frac{\partial}{\partial \bar{r}} (\bar{r}^2 \frac{\partial \bar{\mathcal{E}}}{\partial \bar{r}}) + \bar{r}^2 \bar{\mathcal{B}} (\bar{c} - \bar{\mathcal{C}}) = 0. \tag{8.4}$$

The boundary conditions for the NSC theory are sketched in Fig. 8.1. At the surface, except for the natural boundary condition (4.47), a spatially independent mass flux is chosen, which has been introduced in Equation (5.30). In addition, at the particle center, the boundary conditions

$$\frac{\partial c}{\partial r}(0,t) = 0,\tag{8.5}$$

$$\frac{\partial \hat{c}}{\partial r}(0, t) = 0\tag{8.6}$$

are imposed. Here, Equations (8.5) and (8.6) are needed to ensure the spherical symmetry. What is more, these two boundary conditions ensure the physical requirement that the flux at the particle center vanishes (see Appendix A.4). Finally, the initial conditions are given by

$$c\left(r,0\right) = 0,\tag{8.7}$$

$$
\hat{c}\left(r,0\right) = 0.\tag{8.8}
$$

The cathode material LMO is chosen for the simulations. We consider a typical value of *R*<sup>0</sup> = 1µ*m* for the particle radius. The material parameters for the cathode material LMO have been summarized in Table 7.1, and the parameters ¯β and *l* will be discussed in the following section.

### **8.2 Interface evolution**

Here, we will consider quasistatic insertion of lithium into a particle of LMO at *C* = 0.001. First, we consider six different cases of ¯β and *l* with the same value of the product λ = ¯β*l* <sup>2</sup> = 7×10−<sup>18</sup> *m* 2 . Figs. 8.2 and 8.3 show the normalized lithium concentration ¯*c* and the normalized nonlocal concentration ¯*c*ˆ along the radial coordinate at a state of average concentration of *cavg* = 0.5 for different cases of ¯β and *l*. In Fig. 8.2 , we find that the lithium concentration plots of different cases are almost the same except for the first case ( ¯β = 0.001 and *l* = 83.7 *nm*) which shows a homogeneous state. All other cases exhibit an extremely small difference in the interface situated between the two phases as shown in the inset. Also, Fig. 8.3 shows almost the same nonlocal concentration plots for the different cases. As illustrated in Fig. 8.4, for the last case ( ¯β = 1000 and *l* = 0.084 *nm*), the lithium concentration is the same as the nonlocal concentration, both of which are also identical to the lithium concentration from the Cahn-Hilliard model. According to Eq. (5.46), for the limit case *l* −→ 0 the nonlocal concentration is equal to the lithium concentration ¯*c*, meaning that the Cahn-Hilliard model is obtained and the full nonlocal effect is reduced to a weak one. Therefore, ¯β = 1000 in our simulations is enough to approach the results from the Cahn-Hilliard theory. This is also confirmed by Di Leo et al. [94].

Figure 8.2: Normalized lithium concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 during insertion for different cases of ¯β and *l* [98].

Figure 8.3: Normalized nonlocal concentration ¯*c*ˆ versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 during insertion for different cases of ¯β and *l* [98].

Figure 8.4: Comparison of the concentration plots at *cavg* = 0.5 during insertion from the NSC theory and the Cahn-Hilliard theory, where ¯β = 1000, *l* = 0.084 *nm*, and λ = 7×10−<sup>18</sup> *m* 2 . [98]

Now we focus on the interface evolution shown in the insets of Figs. 8.2 and 8.3. As *l* increases, the interface of the lithium concentration field becomes steeper, but the nonlocal concentration field gets more diffuse, which means that the spatial interaction volume which contributes to the nonlocal effect is larger. The reason is that due to the fixed product λ, the increase of the spatial interaction volume is compensated by a shrinking of the interface which is governed by the characteristic interface length scale *l* [113]. Therefore, the nonlocality makes the interface of the local concentration field ¯*c* steeper. In addition, we should notice that the first case ( ¯β = 0.001 and *l* = 83.7 *nm*) shows no phase segregation, which is explained by Fig. 8.5. Here, zero system penalty energy corresponds to a homogeneous system. For the fixed product λ, as *l* increases, the penalty energy increases. When *l* reaches 83.7 *nm*, the cost of the penalty energy Ψ¯ *penalty* at a phase segregated state would be too high so that as a result of the competition between the different contributions to the total system energy the system stays in the homogenous state in order to minimize the total system free energy. Therefore, the penalty energy vanishes in this case, and it can be seen that the contribution of the penalty energy to the total system free energy can lead to suppression of phase segregation.

Figure 8.5: Normalized penalty energy Ψ¯ *penalty* versus average concentration *cavg* during insertion for different cases of ¯β and *l* [98].

Next, we discuss the evolution of the interface thickness *d* of the lithium concentration field, the interface thickness *dnon* of the nonlocal concentration field, and the characteristic interface length scale *l* with increasing ¯β during insertion, as shown in Fig. 8.6. Here, *d* or *dnon* is defined as the distance over which the concentration deviates 0.05 from the two binodal concentrations of the lithium concentration field or the nonlocal concentration field [113, 141]. We find that *d* increases with increasing ¯β, while both *dnon* and *l* decrease. As ¯β increases, *d* is more close to *dnon*, and they are the same when ¯β reaches 1000 as shown in the inset, when approaching the solution from the Cahn-Hilliard theory as discussed before. However, the difference |*d* −*l*| reduces as ¯β decreases. When ¯β is 1.5, *l* and *d* are equal to 2.16 *nm* and 5.66 *nm*, respectively, and thus *l* is on the order of magnitude of *d*. Therefore, in view of the interpretation of *l* as the characteristic interface length scale, we can regard the two values 1.5 and 2.16 *nm* as reasonable rough estimates of ¯β and *l* for the cathode material LMO, respectively, when fixing λ = 7×10−<sup>18</sup> *m* 2 , which leads to a lower interface thickness of 5.66 *nm* in the NSC theory compared to the interface thickness of 8 *nm* in the Cahn-Hilliard theory.

Figure 8.6: Evolution of the interface thickness *d* of the lithium concentration field, the interface thickness *dnon* of the nonlocal concentration field, and the characteristic interface length scale *l* at *cavg* = 0.5 according to the NSC theory with λ = 7×10−<sup>18</sup> *m* <sup>2</sup> during insertion [98].

## **8.3 Comparison with Cahn-Hilliard model**

As the truly nonlocal NSC theory is an extension of the Cahn-Hilliard theory, we will now compare the results from the two theories. First, we consider the electrode particle as being microsized, say, *R*<sup>0</sup> = 1 µ*m*. Then, as illustrated in Fig. 8.7, the lithium concentration plots obtained from the two theories are almost the same. This means that, for a microsized electrode particle, the nonlocal effect is very weak, and the Cahn-Hilliard theory yields the almost same results as the NSC theory. Therefore, the Cahn-Hilliard theory is well capable of describing the motion of the sharp interface between the two phases for a microsized electrode particle. However, when the electrode particle radius is considered as being nanosized, say, *R*<sup>0</sup> = 100 *nm*, we find in Fig. 8.8 that the lithium concentration at the interface between the two phases obtained from the Cahn-Hilliard theory is different from that of the NSC theory. The interface of the lithium concentration field ¯*c* considering the truly nonlocal effect is steeper than that from the Cahn-Hilliard theory, and the nonlocal concentration field ¯*c*ˆ is the most diffuse among the three concentration plots. Indeed, the ratio of the characteristic interface length scale to the particle radius *l*/*R*<sup>0</sup> in the nanoparticle is larger than that in the microparticle, and, as a result, the larger ratio of the spatial interaction volume to the particle volume contributes to a stronger nonlocal effect.

Figure 8.7: Comparison of the concentration plots at *cavg* = 0.5 during insertion from the NSC theory and the Cahn-Hilliard theory, where *R*<sup>0</sup> = 1 µ*m*, ¯β = 1.5, *l* = 2.16 *nm*, and λ = 7×10−<sup>18</sup> *m* 2 [98].

Figure 8.8: Comparison of the concentration plots at *cavg* = 0.5 during insertion from the NSC theory and the Cahn-Hilliard theory, where *R*<sup>0</sup> = 100 *nm*, ¯β = 1.5, *l* = 2.16 *nm*, and λ = 7×10−<sup>18</sup> *m* 2 [98].

Next, we consider a microparticle with a large microstructure scale. To this end, we choose λ = 7 × 10−<sup>15</sup> *m* 2 , to vary for the sake of comparison the interface thickness in both, NSC and Cahn-Hilliard theory. For fixed ¯β = 1.5, we obtain *l* = 68.31 *nm*. Fig. 8.9 shows the results from the two theories. It is found that the lithium concentration from the Cahn-Hilliard theory is clearly different from that of the NSC theory. The lithium concentration field from the Cahn-Hilliard theory is more diffuse than that from the NSC theory, and the interface of the nonlocal concentration field is very thick, which means that the spatial interaction volume is very large. This can be explained by the increase of the characteristic interface length scale *l* enlarging the spatial interaction volume. As a result, the nonlocal effect is very strong.

Figure 8.9: Comparison of the concentration plots at *cavg* = 0.5 during insertion from the NSC theory and the Cahn-Hilliard theory, where *R*<sup>0</sup> = 1 µ*m*, ¯β = 1.5, *l* = 68.31 *nm*, and λ = 7×10−<sup>15</sup> *m* 2 [98].

## **8.4 Influence of the penalty energy coefficient** β

As the interface thickness *d* is determined by the two parameters ¯β and *l* in the NSC theory, we will study their influence. First, we vary the value of ¯β but fix *l* (*l* = 2.16 *nm*). Figs. 8.10 and 8.11 show the effect of ¯β on the normalized lithium concentration ¯*c* and the normalized nonlocal concentration ¯*c*ˆ at a state of average concentration of *cavg* = 0.5. We find that both, *d* and *dnon*, increase as ¯β increases, but a homogeneous state is obtained for the last case ( ¯β = 15000). This can be explained with the help of Figs. 8.12 and 8.13. Fig. 8.12 shows that with increasing values of ¯β, the interface thicknesses *d* and *dnon* of ¯*c* and ¯*c*ˆ, respectively, tend towards each other asymptotically. This in turn leads to the plots of ¯*c* and ¯*c*ˆ to coincide more and more and, by the way, to approaching the solution of the Cahn-Hilliard theory. The results show that the term *c*¯− ¯*c*ˆ 2 dominates the trend of increasing ¯β in the penalty energy density (4.51). As a result, the system penalty energy in Fig. 8.13 is gradually reduced with increasing value of ¯β. Finally, when ¯β reaches 15000, the system penalty energy vanishes, and a homogeneous state occurs. As discussed before, the asymptotic coincidence of the NSC theory with the Cahn-Hilliard theory for large values of ¯β means a likewise asymptotic reduction to weak nonlocality.

Figure 8.10: Normalized lithium concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 during insertion for different ¯β, where *l* = 2.16 *nm* [98].

Figure 8.11: Normalized nonlocal concentration ¯*c*ˆ versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 during insertion for different ¯β, where *l* = 2.16 *nm* [98].

Figure 8.12: Evolution of the interface thickness *d* of the lithium concentration field, the interface thickness *dnon* of the nonlocal concentration field and the characteristic interface length scale *l* at *cavg* = 0.5 during insertion, where *l* = 2.16 *nm* [98].

Figure 8.13: Normalized penalty energy Ψ¯ *penalty* versus average concentration *cavg* during insertion for different ¯β, where *l* = 2.16 *nm* [98].

## **8.5 Influence of the characteristic interface length scale** *l*

Next, we vary the value of *l* but fix ¯β = 1.5. Figs. 8.14 and 8.15 show the effect of *l* on the normalized lithium concentration ¯*c* and the normalized nonlocal concentration ¯*c*ˆ at a state of average concentration of *cavg* = 0.5. As it has to be expected, increasing the characteristic interface length scale *l* broadens *d* and *dnon*, and increases the difference between ¯*c* and ¯*c*ˆ, which indicates that the nonlocal effect is enhanced. This is consistent with the increase of the maxima of the system penalty energy in Fig. 8.16 as *l* increases. Interestingly, comparing the case ( ¯β = 1.5, λ = 7×10−<sup>14</sup> *m* 2 ) in Fig. 8.14 to the case (*l* = 2.16 *nm*, λ = 7×10−<sup>14</sup> *m* 2 ) in Fig. 8.10, we find that it is a phase segregated state for the former case but a homogeneous state for the latter case, although the product λ is the same. Indeed, ¯β in the latter case is far more than 1000, which means that the result from the Cahn-Hilliard theory is approached. On the other hand, the former case represents a result reasonably accounting for a truly nonlocal effect due to the fact that the characteristic interface length scale *l* is on the order of magnitude of the interface thickness *d*. Therefore we conclude that, the Cahn-Hilliard theory is not capable of describing the interface evolution of a microstructure possessing a length scale comparable to the total system size. Now, we focus on the evolution of *d* and *dnon* with increasing *l* for fixed ¯β, as shown in Fig. 8.17. It is found that the difference |*d* −*dnon*| initially increases as *l* increases. Beyond *l* = 68.3 *nm*, the difference |*d* −*dnon*| gets smaller. Finally, for the case of *l* = 683 *nm* and λ = 7×10−<sup>13</sup> *m* 2 , which is not shown in the figures, there is a homogeneous state and, as a consequence both *d* and *dnon* vanish with the trivial consequence |*d*−*dnon*| = 0. The reason for this change of the trend of |*d* − *dnon*| lies in the competition of the different terms contributing to the system free energy. As the system penalty energy is about to get too large (see Fig. 8.16), a homogenous solution yields the minimum system free energy. In other words, too large a penalty energy suppresses phase segregation. In addition, *d* is always on the order of magnitude of *l* in the different cases. Therefore, the value 1.5 can be confirmed again as a reasonable rough estimate of ¯β for the NSC theory.

Figure 8.14: Normalized lithium concentration ¯*c* versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 during insertion for different *l*, where ¯β = 1.5 [98].

Figure 8.15: Normalized nonlocal concentration ¯*c*ˆ versus normalized radial coordinate *r*/*R*<sup>0</sup> at *cavg* = 0.5 during insertion for different *l*, where ¯β = 1.5 [98].

Figure 8.16: Normalized penalty energy Ψ¯ *penalty* versus average concentration *cavg* during insertion for different *l*, where ¯β = 1.5 [98].

Figure 8.17: Evolution of the interface thickness *d* of the lithium concentration field, the interface thickness *dnon* of the nonlocal concentration field and the interfacial characteristic interface length scale *l* at *cavg* = 0.5 during insertion, where ¯β = 1.5 [98].

# **9 Conclusion**

The occurence of phase segregation within electrode particles, which is considered as a resource of large mechanical stresses, can be a possible origin for the battery capacity fade. Thus, it is the objective of this work to study phase changes and mechanical stresses in electrode particles by means of phasefield simulations. This work was first dedicated to the understanding of phase changes, mechanical deformation, and stress evolution in NaFPO electrode particles of NIBs. Furthermore, the particle size and average concentration dependent miscibility gap for the three cathode materials LMO, LiFPO, and NaPPO were studied. Finally, a NSC theory, which is an extension of the Cahn-Hilliard theory, was introduced.

We formulated a non-linear Cahn-Hilliard type model with mechanics. For the mechanical part, besides SDT two different finite deformation elasticity formulations were introduced and compared. This work has addressed the theoretical formulation, numerical implementation, and application of the mechanically coupled phase-field theory.

To start with, a phase-field model for the cathode material NaFPO of NIBs was studied for the first time, using a phase-field theory coupling the Cahn-Hilliard equation to finite deformation elasticity. As a major novelty, the material parameters for NaFPO were determined. For example, α1, α2, and λ, all of which are the key parameters in the phase-field model, were determined. The determination of these key parameters provides a significant input for the future phase-field work for NaFPO. We implemented the fourth-order nonlinear initial-boundary-value problem of the model in COMSOL Multiphysics<sup>r</sup> to solve by the finite element method the spherically symmetric problem of sodium insertion into or extraction from a NaFPO particle of NIBs. Our model captures the important feature that distinguishes NaFPO from LiFPO, i.e, phase segregation into a sodium-poor phase *FePO*<sup>4</sup> and a sodium-rich phase *Na*2/3*FePO*4. The difference for the concentration and stress plots between SDT and the finite deformation theories can not be neglected for NaFPO. In particular, the stresses in the two phases are higher for SDT. The above difference suggests that the finite deformation elasticity is preferred for the cathode material NaFPO. This is different from that SDT has a sufficient capacity to represent the deformation of LiFPO. Although the results from the two finite deformation theories are almost the same, according to our experience, elasticity based on logarithmic elastic strain is numerically more efficient than elasticity based on elastic Green strain. On the other hand, we found that, for the same stage of average concentration, the miscibility gap of NaFPO is smaller than that of LiFPO, but the stresses in a NaFPO particle during phase segregation are larger, which may explain the existence of a wide range of solid solution *NaxFePO*<sup>4</sup> for 2/3 < *x* < 1 to avoid even higher stresses in the NaFPO particle. In addition, the suppression of phase segregation by the elastic strain energy is more pronounced in NaFPO compared to LiFPO. Here, we just focus on phase changes in the two-phase region of NaFPO (0 < *x* < 2/3), ignoring the single-phase region (2/3 < *x* < 1). As an outlook, the current model may be extended into the whole region (0 < *x* < 1). As both the partial molar volume and Young's modulus are regarded as constants here, another major subject of future study is the effect of the concentration dependence of these two quantities for NaFPO.

For the second aim of this work, by means of a phase-field theory coupling the Cahn-Hilliard equation to finite deformation elasticity, we studied the dependence of the miscibility gap on particle size and average concentration for the three cathode materials LMO, LiFPO, and NaFPO. We found that, a reduced miscibility gap in the center of the spinodal region of average concentration is found for radii below 55.5 *nm*, 20.5 *nm*, and 17.5 *nm* for the three cathode materials LMO, LiFPO, and NaFPO, respectively. The miscibility gap shrinks to accommodate in particles of reduced size the gradient of concentration in the interface between phases. The critical radius below which phase segregation is inhibited for LMO is *R<sup>c</sup>* = 10.9 *nm*, which is larger than the critical radii of 6 *nm* for LiFPO and 6.4 *nm* for NaFPO. Concerning the evolution of the miscibility gap during the insertion process, it was found that for smaller particles it increases as average concentration increases. However, the miscibility gap is constant during the whole process of insertion for radii above 125 *nm*, 73 *nm*, and 75 *nm* for the three cathode materials LMO, LiFPO, and NaFPO, respectively. The average concentration dependent miscibility gap is physically explained by the minimization of system free energy. On the other hand, when mechanics is accounted for, minimization of the total energy leads to a reduction of the miscibility gap and a widened interface region. Among the three investigated cathode materials, mechanics has the strongest influence on LMO on the miscibility gap, while the tensile stresses in this material are the smallest. The elastic strain energy leads to a larger threshold value of the particle radius for the size independent miscibility gap, and a larger critical particle size below which phase segregation is completely inhibited. Our work suggests that one way to improve the mechanical stability of cathode materials is that the miscibility gap is suppressed by choosing sufficiently small particle sizes according to the respective material. For the future challenge to design nanoscale insertion materials, LMO possesses the largest particle size for expanding the solid solution range in which particle stresses are eliminated among the three investigated cathode materials. Here, neither surface tension nor surface wetting [50, 51, 142] is taken into account. According to Stein et al. [142], surface tension plays an increasingly unnegligible role in the electro-chemo-mechanical behavior of electrode nanoparticles as the particle size decreases. Investigating the effect of surface tension on the miscibility gap in phase separating insertion materials can be subject of future work. In order to avoid computationally expensive three-dimensional simulations, the cathodic particles are assumed to be of spherical symmetry under spherically symmetric boundary conditions. However, the particle shape has an important influence on the phase segregation and the stress state [88, 95, 116]. According to Zhao et al. [95], the core-shell phase segregation may not occur due to the geometric anisotropy in plate- and needle-like particles, and the stress magnitudes in the particles of the above two types can be different compared to those in particles of spherical symmetry. Investigating the effect of the particle shape may also be subject of future work.

Finally, a NSC theory was introduced from a general formulation of a nonlocal free energy density by using the Green's function as weight function to describe diffusion and phase changes in a material showing phase segregation. This theory was shown to be a generalization of the classical Cahn-Hilliard theory. In contrast to the Cahn-Hilliard theory which is weakly nonlocal, the NSC theory is truly nonlocal and is based not only on the species concentration *c* but also on the nonlocal species concentration ˆ*c*. The NSC theory is governed by two coupled second-order PDEs which are amenable to being solved with a standard finite element implementation and, thus, is computationally less demanding than the fourth-order Cahn-Hilliard equation. In the NSC theory, the nonlocal free energy density is split into two parts, which are the penalty energy density ψ *penalty* and the variance energy density ψ *variance* . The interface thickness *d* is controlled by two independent parameters, namely the penalty energy coefficient ¯β and the characteristic interface length scale *l*. We implemented the NSC theory in COMSOL Multiphysics<sup>r</sup> to solve by the finite element method the spherically symmetric problem of lithium insertion into a LMO particle of a lithium ion battery. It was found that ¯β = 1.5 and *l* = 2.16 *nm* can be regarded as reasonable rough estimates for the cathode material LMO. The nonlocality makes the interface of the local concentration field steeper, when compared to the Cahn-Hilliard theory. Also, the NSC theory is more accurate than the Cahn-Hilliard theory for a nanoparticle of an electrode material or, more general, a comparable scale of the microstructure in relation to the system size. A high system penalty energy can lead to suppression of phase segregation. On the other hand, both interface thicknesses *d* and *dnon*

increase as ¯β or *l* increase. Furthermore, increasing ¯β weakens the nonlocal effect while increasing *l* enhances the nonlocality.

# **Bibliography**


domino-cascade model," *Nature Materials*, vol. 7, no. 8, pp. 665–671, 2008.


"Accommodating high transformation strains in battery electrodes via the formation of nanoscale intermediate phases: operando investigation of olivine *NaFePO*4," *Nano Letters*, vol. 17, no. 3, pp. 1696–1702, 2017.


deformations: application to phase-separating Li-ion electrode materials," *Journal of the Mechanics and Physics of Solids*, vol. 70, pp. 1–29, 2014.


behavior," *Continuum Mechanics and Thermodynamics*, vol. 26, no. 4, pp. 551–562, 2014.


[142] P. Stein, Y. Zhao, and B.-X. Xu, "Effects of surface tension and electrochemical reactions in li-ion battery electrode nanoparticles," *Journal of Power Sources*, vol. 332, pp. 154–169, 2016.

# **List of Figures**









# **List of Tables**


# **A Appendix**

## **A.1 Helmholtz equation (4.46) as governing partial differential equation for nonlocal species concentration**

The Green's function is the solution of the Helmholtz equation (4.46) with the source term replaced by the Dirac function δ(~*x*−~*y*), i.e., the solution of

$$G(\vec{x}; \vec{\mathfrak{y}}) - l^2 \nabla^2 G(\vec{x}; \vec{\mathfrak{y}}) = \mathcal{S}(\vec{x} - \vec{\mathfrak{y}}),\tag{A.1.1}$$

for the boundary condition (4.47). Integrating Equation (A.1.1) over the problem domain B with respect to spatial position~*y*, one obtains

$$\begin{split} \int\_{\beta \overline{\theta}} G(\vec{x}; \vec{y})dV &= \int\_{\beta \overline{\theta}} \delta(\vec{x} - \vec{y})dV + l^2 \int\_{\beta \overline{\theta}} \nabla^2 G(\vec{x}; \vec{y})dV \\ &= \quad 1 + l^2 \int\_{\partial \beta \overline{\theta}} \vec{\nabla} G(\vec{x}; \vec{y}) \cdot \vec{n} dA \\ &= \quad 1, \end{split} \tag{A.1.2}$$

where the term related to *l* <sup>2</sup> vanishes by the divergence theorem and the boundary condition (4.47). Consequently, if the Greens's function *G*(~*x*;~*y*) is continued by zero values outside the problem domain B, it satisfies the natural weight function property (4.28), and, thus it is an admissible choice of the weight function, i.e. ω(~*y*;~*x*) = *G*(~*y*;~*x*). Notice that *G*(~*y*;~*x*) = *G*(~*x*;~*y*) since Green's functions are symmetric for linear problems with constant coefficients [114]. Thus, we introduce the nonlocal species concentration according to Equation (4.35) as

$$
\tilde{c}(\vec{x}) \quad = \frac{A}{c\_{\text{max}} \mathfrak{B}} \int\_{\mathcal{A}} \tilde{c}(\vec{y}) G(\vec{y}; \vec{x}) dV. \tag{A.1.3}
$$

In view of the choice of *G* as the weight function, the definition (4.36) together with Equation (A.1.2) yields

$$\mathcal{B} = \frac{A}{c\_{\text{max}}} \int\_{\mathcal{B}} G(\vec{\mathbf{y}}; \vec{\mathbf{x}}) dV = \frac{A}{c\_{\text{max}}},\tag{A.1.4}$$

such that Equation (A.1.3) simplifies to

$$
\tilde{\mathfrak{E}}(\vec{\mathfrak{x}}) \quad = \int\_{\beta \vec{\mathcal{B}}} \tilde{c}(\vec{\mathfrak{y}}) G(\vec{\mathfrak{y}}; \vec{\mathfrak{x}}) dV. \tag{A.1.5}
$$

By the definition of the Dirac function

$$
\vec{c}(\vec{x}) \quad = \int\_{\beta\overline{\theta}} \mathcal{S}(\vec{x} - \vec{y}) \vec{c}(\vec{y}) dV \tag{A.1.6}
$$

holds. Now, substituting Equation (A.1.1) into Equation (A.1.6) gives

$$\begin{split} \bar{c}(\vec{\boldsymbol{x}}) &= \int\_{\mathcal{A}} \bar{c}(\vec{\boldsymbol{y}}) G(\vec{\boldsymbol{x}}; \vec{\boldsymbol{y}}) dV - \int\_{\mathcal{A}} l^{2} \bar{c}(\vec{\boldsymbol{y}}) \nabla^{2} G(\vec{\boldsymbol{x}}; \vec{\boldsymbol{y}}) dV \\ &= \int\_{\mathcal{A}} \bar{c}(\vec{\boldsymbol{y}}) G(\vec{\boldsymbol{y}}; \vec{\boldsymbol{x}}) dV - l^{2} \nabla^{2} \int\_{\mathcal{A}} \bar{c}(\vec{\boldsymbol{y}}) G(\vec{\boldsymbol{y}}; \vec{\boldsymbol{x}}) dV \\ &= \quad \bar{\mathcal{E}}(\vec{\boldsymbol{x}}) - l^{2} \nabla^{2} \tilde{\mathcal{E}}(\vec{\boldsymbol{x}}). \end{split} \tag{A.1.7}$$

Therefore, the inhomogenous Helmholtz equation (4.46) is the governing equation for nonlocal species concentration when ¯*c*ˆ is defined according to the integral representation (4.35) with the Green's function as the weight function. As an additional implication, as also shown in [114], the coefficients related to derivatives of fourth-order and higher of ¯*c* in Equation (4.45) need to vanish, such that equation (4.45) simplifies to the Helmholtz equation (4.46).

# **A.2 Derivation of the system free energy in the nonlocal species concentration theory**

The free energy density (4.27) in the nonlocal model is reformulated as

$$\begin{split} \Psi(\vec{x}) &= \quad \Psi^{\mu\nu\rho} \left( \vec{c}(\vec{x}) \right) + \frac{1}{2} A \int\_{\beta\theta} \mathfrak{o}(\rho) \left( \vec{c}(\vec{x}) - \vec{\bar{c}}(\vec{x}) \right)^{2} dV \\ &+ \frac{1}{2} A \int\_{\beta\theta} \mathfrak{o}(\rho) \left( \vec{c}(\vec{x}) - \vec{c}(\vec{y}) \right)^{2} - \mathfrak{o}(\rho) \left( \vec{c}(\vec{x}) - \vec{\bar{c}}(\vec{x}) \right)^{2} dV. \end{split} \tag{A.2.1}$$

Note that ~*y* is the variable of spatial integration. Using Equation (4.36), the second term on the right hand side of Equation (A.2.1) gives

$$\begin{split} \frac{1}{2}A \int\_{\left(\vec{\mathcal{B}}\right)} \mathcal{a}(\boldsymbol{\rho}) \left(\vec{c}(\vec{\boldsymbol{x}}) - \vec{\mathcal{E}}(\vec{\boldsymbol{x}})\right)^{2}dV &= \frac{1}{2}c\_{\max} \mathcal{B}\left(\vec{c}(\vec{\boldsymbol{x}}) - \vec{\mathcal{E}}(\vec{\boldsymbol{x}})\right)^{2} \\ &= \quad \Psi^{penalty}(\vec{c}(\vec{\boldsymbol{x}}), \vec{\mathcal{E}}(\vec{\boldsymbol{x}})). \end{split} \tag{A.2.2}$$

Using Equations (4.35) and (A.1.2), the third term on the right hand side of Equation (A.2.1) then yields

$$\frac{1}{2}A\int\_{\beta\overline{\theta}}\mathfrak{o}(\mathfrak{p})\left(\bar{c}(\vec{x})-\bar{c}(\vec{y})\right)^{2}-\mathfrak{o}(\mathfrak{p})\left(\bar{c}(\vec{x})-\bar{\tilde{c}}(\vec{x})\right)^{2}dV$$

$$=(A-c\_{\max}\mathcal{B})\bar{c}(\vec{x})\bar{\tilde{c}}(\vec{x})-\frac{1}{2}A\bar{\tilde{c}}(\vec{x})^{2}+\frac{1}{2}A\int\_{\beta\overline{\theta}}\mathfrak{o}(\mathfrak{p})\bar{c}(\vec{y})^{2}dV.\tag{A.2.3}$$

The first term on the right hand side of Equation (A.2.3) vanishes by Equation (A.1.4), and Equation (A.2.3) then can be manipulated further:

$$\begin{split} & \quad -\frac{1}{2}A\tilde{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{x}})^{2} + \frac{1}{2}A \int\_{\mathcal{\mathcal{B}}} \boldsymbol{\mathfrak{o}}(\rho)\boldsymbol{\bar{\varepsilon}}(\vec{\boldsymbol{y}})^{2}dV \\ &= \quad \quad \quad \frac{1}{2}A \Big( \int\_{\mathcal{\mathcal{B}}} \boldsymbol{\mathfrak{o}}(\rho)\boldsymbol{\bar{\varepsilon}}(\vec{\boldsymbol{y}})^{2}dV - \bar{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{x}})^{2} \Big) \\ &= \quad \quad \quad \frac{1}{2}A \Big( \int\_{\mathcal{\mathcal{B}}} \boldsymbol{\mathfrak{o}}(\rho)\boldsymbol{\bar{\varepsilon}}(\vec{\boldsymbol{y}})^{2}dV + \bar{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{x}})^{2} - 2\bar{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{x}})^{2} \Big) \\ &= \quad \quad \quad \frac{1}{2}A \Big( \int\_{\mathcal{\mathcal{B}}} \boldsymbol{\mathfrak{o}}(\rho)\boldsymbol{\bar{\varepsilon}}(\vec{\boldsymbol{y}})^{2}dV + \bar{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{x}})^{2} - 2\bar{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{x}}) \int\_{\mathcal{\mathcal{B}}} \boldsymbol{\mathfrak{o}}(\rho)\boldsymbol{\bar{\varepsilon}}(\vec{\boldsymbol{y}})dV \Big) \\ &= \quad \quad \quad \quad \frac{1}{2}c\_{\text{max}}\mathcal{B} \int\_{\mathcal{\mathcal{B}}} \boldsymbol{\mathfrak{o}}(\rho)\left(\bar{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{y}}) - \bar{\boldsymbol{\varepsilon}}(\vec{\boldsymbol{x}})\right)^{2}dV \\ &= \quad \quad \quad \quad \quad \$$

As shown in appendix A.1, the Green's function *G*(~*y*;~*x*) of the Helmholtz equation (4.46) satisfies the natural weight function property (4.28). Thus, it is chosen as the weight function of the NSC theory. As a result, we obtain the system free energy (4.50) in which the nonlocal free energy density is composed of two terms, namely the penalty energy density ψ *penalty* and variance energy ψ *variance* .

# **A.3 Derivation of the final coupled diffusion equation of the spherically symmetric boundary value problem**

## **A.3.1 The coupled Cahn-Hilliard equation with elastic Green strain**

Since Gradµ *cp R* (*cR*,E *e* ) is extremely complicated and long, it makes sense to use the following substitutions:

$$\mathbf{s} \quad = \quad \mathfrak{Q}(c\_R - c\_0), \tag{A.3.1}$$

$$q\_{\parallel} = \begin{array}{c} \frac{\partial u\_R}{\partial R}, \\\\ \end{array} \tag{A.3.2}$$

$$
\omega\_{\parallel} = \begin{array}{c c}
\underline{u\_R} \\
\underline{R}
\end{array}
\tag{A.3.3}
$$

$$m\_{\parallel} = \frac{\mathbf{v}}{1 - 2\mathbf{v}},\tag{A.3.4}$$

$$X\_{\
u} = \left(1 + s\right)^{-\frac{8}{3}},\tag{A.3.5}$$

$$Y\_{\alpha} = \begin{array}{c} \left(1+s\right)^{-\frac{5}{3}}, \\\\ \left(1+s\right)^{-\frac{2}{3}} \end{array} \tag{A.3.6}$$

$$Z\_{\;\;\;\;\phi} = \;\;(1+s)^{-\frac{\sharp}{3}},\tag{A.3.7}$$

$$\mathcal{Q}\_{-}=\left(\mathbb{I}+q\right)^{2},\tag{A.3.8}$$

$$W \quad = \quad \left(1+w\right)^2. \tag{A.3.9}$$

According to Table 5.1, all the above introduced parameters (A.3.1 - A.3.9) are dimensionless.

For the gradient of the chemical potential, Gradµ *mwp R* and Gradµ *gd R* have been expressed in Equations (5.25) and (5.26), respectively. Again, here *c* and *r* in Equations (5.25) and (5.26) should be respectively replaced by *c<sup>R</sup>* and *R*. Using Equation (5.80), Gradµ *cp R* is expressed as

$$\begin{split} \text{Grad}\mu\_{\tilde{\mathbf{g}}}^{\mu\prime\prime}(c\_{R},\mathbf{E}^{\prime}) &= \qquad G\Omega \left[\frac{1}{2}(ZQ-1)(-\frac{2}{3}\frac{\partial}{\partial R}QY + 2(1+q)Z\frac{\partial}{\partial R}) \\ &+ (-\frac{2}{3}YW\frac{\partial\overline{s}}{\partial R} + 2Z(\frac{q}{R}-\overline{w})(1+w)(WZ-1)) \\ &+ 2w(-1+\frac{1}{2}(-1+QZ)+WZ) \\ &\cdot (-\frac{1}{3}ZY\frac{\partial}{\partial R} + Z(1+q)\frac{\partial}{\partial R} - \frac{2}{3}\gamma W\frac{\partial}{\partial R} \\ &+ 2Z(1+w)(\frac{q}{R}-\overline{w})\Big{)}\delta\_{\overline{R}} \\ &+ G\Omega(1+w)\left(\frac{\overline{s}}{2}(Q-1)QX\frac{\partial}{\partial R}\right) \\ &- \frac{1}{3}ZY\left(-\frac{2}{3}Y\frac{\partial}{\partial R}Q + 2Z(1+q)\frac{\partial}{\partial R}\right) \\ &- \frac{2}{3}(ZQ-1)(1+q)\frac{\partial}{\partial R}Y + \frac{10}{9}XW\frac{\partial}{\partial R}(WZ-1) \\ &- \frac{4}{3}(WZ-1)Y(1+w)(\frac{q}{R}-\frac{w}{R}) - \frac{2}{3}YY(\frac{-2}{3})YW\frac{\partial}{\partial R} \\ &+ 2(w+1)2(\frac{\overline{\dot{s}}}{R}-\frac{\dot{W}}{R}) \\ &+ 2w(\frac{\overline{s}}{3}Q\frac{\partial}{\partial R}-\frac{2}{3}(1+q)\frac{\partial}{\partial R}Y + \frac{10}{9}W\frac{\partial}{\partial R}Z \\ &- \frac{4}{3}(1+w)Y(\frac{q}{R}-\frac{\$$

(A.3.10) <sup>182</sup>

Substituting Equations (5.25), (5.26), and (A.3.10) into Equation (5.68), we can derive the final dimensionless coupled diffusion equation

0 = *R*¯ 2 ∂ *c*¯*<sup>R</sup>* ∂*t*¯ + ∂ ∂*R*¯ −*R*¯ 2 (1+α2*c*¯*R*(1−*c*¯*R*))<sup>∂</sup> *<sup>c</sup>*¯*<sup>R</sup>* ∂*R*¯ +*R*¯ <sup>2</sup> λ *R* 2 0 *c*¯*R*(1−*c*¯*R*) · ∂ 3 *c*¯*R* ∂*R*¯<sup>3</sup> + 2 *R*¯ ∂ 2 *c*¯*R* ∂*R*¯<sup>2</sup> − 2 *R*¯2 ∂ *c*¯*<sup>R</sup>* ∂*R*¯ −*R*¯ 2 *c*¯*R*(1−*c*¯*R*) · *G*¯Ω¯ 1 2 (*ZQ*−1)(− 2 3 ∂ *s* ∂*R*¯ *QY* +2(1+*q*)*Z* ∂*q* ∂*R*¯ ) + (− 2 3 *YW* ∂ *s* ∂*R*¯ +2*Z*( *q R*¯ − *w R*¯ )(1+*w*))(*W Z* −1) +2*n*(−1+ 1 2 (−1+*QZ*) +*W Z*) ·(− 1 3 *QY* ∂ *s* ∂*R*¯ +*Z*(1+*q*) ∂*q* ∂*R*¯ − 2 3 *WY* ∂ *s* ∂*R*¯ <sup>+</sup>2*Z*(1+*w*)( *<sup>q</sup> R*¯ − *w R*¯ )) +*G*¯Ω¯ (1+*s*) 5 9 (*ZQ*−1)*QX* ∂ *s* ∂*R*¯ − 1 3 *QY*(− 2 3 *Y* ∂ *s* ∂*R*¯ *Q* +2*Z*(1+*q*) ∂*q* ∂*R*¯ )− 2 3 (*ZQ*−1)(1+*q*) ∂*q* ∂*R*¯ *Y* + 10 9 *XW* ∂ *s* ∂*R*¯ (*W Z* −1) − 4 3 (*W Z* <sup>−</sup>1)*Y*(1+*w*)( *<sup>q</sup> R*¯ − *w R*¯ )− 2 3 *WY*(− 2 3 *YW* ∂ *s* ∂*R*¯ +2(*w*+1)*Z*( *q R*¯ − *w R*¯ )) +2*n*( 5 9 *QX* ∂ *s* ∂*R*¯ − 2 3 (1+*q*) ∂*q* ∂*R*¯ *Y* + 10*W X* 9 ∂ *s* ∂*R*¯ − 4 3 (1+*w*)*Y*( *q R*¯ − *w R*¯ ))(*QZ* <sup>−</sup><sup>1</sup> 2 +*W Z* −1)−2*n*( *QY* 3 + 2*WY* 3 ) ·(− 1 3 *QY* ∂ *s* ∂*R*¯ + (1+*q*) ∂*q* ∂*R*¯ *Z* − 2 3 ∂ *s* ∂*R*¯ *YW* <sup>+</sup>2(*w*+1)( *<sup>q</sup> R*¯ − *w R*¯ )*Z*) +*G*¯ ∂ *s* ∂*R*¯ − Ω¯ 3 (*ZQ*−1)*QY* − 2Ω¯ 3 (*W Z* −1)*WY* +2*n*(− Ω¯ 3 *QY* − 2Ω¯ 3 *WY*)(−1+ 1 2 (−1+*QZ*) +*W Z*) . (A.3.11)

## **A.3.2 The coupled Cahn-Hilliard equation with logarithmic elastic strain**

Since Gradµ *cp R* (*cR*,E *e log*) is extremely complicated and long as well, it makes sense to use the following substitutions:

$$X^\* \quad = \quad (1+s)^{-\frac{7}{3}},\tag{A.3.12}$$

$$Y^\* \quad = \quad (1+s)^{-\frac{4}{3}},\tag{A.3.13}$$

$$Z^\* \quad = \quad \text{(1+s)}^{-\frac{1}{3}}.\tag{A.3.14}$$

According to Table 5.1, all the above introduced parameters (A.3.12 - A.3.14) are dimensionless.

Using Equation (5.92), Gradµ *cp R* (*cR*,E *e log*) is expressed as

Gradµ *cp R* (*cR*,E *e log*) = *G*Ω 2lnλ *e* 1 λ *e* 1 (− 1 3 ∂ *s* ∂*R Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R* + 4lnλ *e* 2 λ *e* 2 (− 1 3 ∂ *s* ∂*R Y* ∗ (1+*w*) +*Z* ∗ ( *q R* − *w R* ) +2*n* lnλ *e* 1 +2lnλ *e* 2 · 1 λ *e* 1 (− 1 3λ *e* 1 ∂ *s* ∂*R Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R* + 2 λ *e* 2 (− 1 3 ∂ *s* ∂*R Y* ∗ (1+*q*) +*Z* ∗ ( *q R* − *w R* ) ~*e<sup>R</sup>* +*G*(1+*s*) " − 2Ω 3 1 λ *e* 1 2 (1+*q*)*Y* ∗ − 1+*q* 3 ∂ *s* ∂*R Y* ∗ +*Z* ∗ ∂*q* ∂*R* − lnλ *e* 1 λ *e* 1 2 (1+*q*)*Y* ∗ − 1 3 ∂ *s* ∂*R Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R* + lnλ *e* 1 λ *e* 1 *Y* ∗ ∂*q* ∂*R* − 4lnλ *e* 1 3λ *e* 1 (1+*q*)*X* ∗ ∂*q* ∂*R* − 4Ω 3 1 λ *e* 2 2 (1+*w*)*Y* ∗ · − 1 3 ∂ *s* ∂*R Y* ∗ (1+*w*) +*Z* ∗ ( *q R* − *w R* ) − lnλ *e* 2 λ *e* 2 2 (1+*w*)*Y* ∗ − *Y* ∗ 3 ∂ *s* ∂*R* (1+*w*) +*Z* ∗ ( *q R* − *w R* ) + lnλ *e* 2 λ *e* 2 ( *q R* − *w R* )*Y* <sup>∗</sup> − 4lnλ *e* 2 3λ *e* 2 (1+*w*)*X* ∗ ∂ *s* ∂*R* 

+2*n* 1 λ *e* 1 − 1 3 ∂ *s* ∂*R Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R* + 2 λ *e* 2 − 1 3 ∂ *s* ∂*R Y* ∗ (1+*q*) +*Z* ∗ ( *q R* − *w R* ) · − Ω 3λ *e* 1 (1+*q*)*Y* <sup>∗</sup> − 2Ω 3λ *e* 2 (1+*w*)*Y* ∗ +2*n* lnλ *e* 1 +2lnλ *e* 2 Ω 3λ *e* 1 2 (1+*q*)*Y* ∗ · − 1 3 ∂ *s* ∂*R Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R* − Ω 3λ *e* 1 ∂*q* ∂*R Y* <sup>∗</sup> + 4Ω 9λ *e* 1 (1+*q*)*X* ∗ ∂ *s* ∂*R* + 2Ω 3λ *e* 2 2 (1+*w*)*Y* ∗ · − *Y* ∗ 3 ∂ *s* ∂*R* (1+*w*) +*Z* ∗ ( *q R* − *w R* ) − 2Ω 3λ *e* 2 ( *q R* − *w R* )*Y* <sup>∗</sup> + 8Ω 9λ *e* 2 (1+*w*)*X* ∗ ∂ *s* ∂*R* ~*e<sup>R</sup>* +*G* ∂ *s* ∂*R* − 2Ω 3 lnλ *e* 1 λ *e* 1 *Y* ∗ (1+*q*) − 4Ω 3 lnλ *e* 2 λ *e* 2 (1+*w*)*Y* <sup>∗</sup> +2*n*Ω lnλ *e* 1 +2lnλ *e* 2 · − 1 3 1 λ *e* 1 *Y* ∗ (1+*q*)− 2 3 1 λ *e* 2 *Y* ∗ (1+*w*) ~*eR*. (A.3.15)

Substituting Equations (5.25), (5.26), and (A.3.15) into Equation (5.68), we can derive the final dimensionless coupled diffusion equation

0 = *R*¯ 2 ∂ *c*¯*<sup>R</sup>* ∂*t*¯ + ∂ ∂*R*¯ −*R*¯ 2 (1+α2*c*¯*R*(1−*c*¯*R*))<sup>∂</sup> *<sup>c</sup>*¯*<sup>R</sup>* ∂*R*¯ +*R*¯ <sup>2</sup> λ *R* 2 0 *c*¯*R*(1−*c*¯*R*) · ∂ 3 *c*¯*R* ∂*R*¯<sup>3</sup> + 2 *R*¯ ∂ 2 *c*¯*R* ∂*R*¯<sup>2</sup> − 2 *R*¯2 ∂ *c*¯*<sup>R</sup>* ∂*R*¯ −*R*¯ 2 *c*¯*R*(1−*c*¯*R*) · *G*¯Ω¯ 2lnλ *e* 1 λ *e* 1 (− 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R*¯ + 4lnλ *e* 2 λ *e* 2 (− 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*w*) +*Z* ∗ ( *q R*¯ − *w R*¯ ) +2*n* lnλ *e* 1 +2lnλ *e* 2 1 λ *e* 1 (− 1 3λ *e* 1 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R*¯ + 2 λ *e* 2 (− 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ( *q R*¯ − *w R*¯ ) +*G*¯(1+*s*) " − 2Ω¯ 3 1 λ *e* 1 2 (1+*q*)*Y* ∗ − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R*¯ − lnλ *e* 1 λ *e* 1 2 (1+*q*)*Y* ∗ − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R*¯ + lnλ *e* 1 λ *e* 1 *Y* ∗ ∂*q* ∂*R*¯ − 4lnλ *e* 1 3λ *e* 1 (1+*q*)*X* ∗ ∂*q* ∂*R*¯ − 4Ω¯ 3 1 λ *e* 2 2 (1+*w*)*Y* ∗ − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*w*) +*Z* ∗ ( *q R*¯ − *w R*¯ ) − lnλ *e* 2 λ *e* 2 2 (1+*w*)*Y* ∗ − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*w*) +*Z* ∗ ( *q R*¯ − *w R*¯ ) + lnλ *e* 2 λ *e* 2 ( *q R*¯ − *w R*¯ )*Y* <sup>∗</sup> − 4lnλ *e* 2 3λ *e* 2 (1+*w*)*X* ∗ ∂ *s* ∂*R*¯ +2*n* 1 λ *e* 1 − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R*¯ 

+ 2 λ *e* 2 − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ( *q R*¯ − *w R*¯ ) · − Ω¯ 3λ *e* 1 (1+*q*)*Y* <sup>∗</sup> − 2Ω¯ 3λ *e* 2 (1+*w*)*Y* ∗ +2*n* lnλ *e* 1 +2lnλ *e* 2 Ω¯ 3λ *e* 1 2 (1+*q*)*Y* ∗ − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*q*) +*Z* ∗ ∂*q* ∂*R*¯ − Ω¯ 3λ *e* 1 ∂*q* ∂*R*¯ *Y* <sup>∗</sup> + 4Ω¯ 9λ *e* 1 (1+*q*)*X* ∗ ∂ *s* ∂*R*¯ + 2Ω¯ 3λ *e* 2 2 (1+*w*)*Y* ∗ − 1 3 ∂ *s* ∂*R*¯ *Y* ∗ (1+*w*) +*Z* ∗ ( *q R*¯ − *w R*¯ ) − 2Ω¯ 3λ *e* 2 ( *q R*¯ − *w R*¯ )*Y* <sup>∗</sup> + 8Ω¯ 9λ *e* 2 (1+*w*)*X* ∗ ∂ *s* ∂*R*¯ +*G*¯ ∂ *s* ∂*R*¯ − 2Ω¯ 3 lnλ *e* 1 λ *e* 1 *Y* ∗ (1+*q*)− 4Ω¯ 3 lnλ *e* 2 λ *e* 2 (1+*w*)*Y* ∗ 2*n*Ω¯ lnλ *e* 1 +2lnλ *e* 2 · − 1 3 1 λ *e* 1 *Y* ∗ (1+*q*)− 2 3 1 λ *e* 2 *Y* ∗ (1+*w*) . (A.3.16)

# **A.4 Lithium flux at the particle center for the spherically symmetric boundary value problem of the nonlocal species concentration theory**

Here, we will demonstrate that the two boundary conditions (8.5) and (8.6) ensure that the flux at the particle center vanishes. This is a physical requirement in order to guarantee the conservation of lithium matter within the particle.

Using Equations (4.4), (4.37), and (4.10), we can obtain the normalized lithium flux in the one-dimensional particle model as

$$\vec{J} = -\left(1 + \mathcal{Q}\_2 \bar{c} (1 - \bar{c})\right) \frac{\partial \bar{c}}{\partial \bar{r}} - \bar{\beta} \bar{c} (1 - \bar{c}) (\frac{\partial \bar{c}}{\partial \bar{r}} - \frac{\partial \bar{\hat{c}}}{\partial \bar{r}}),\tag{A.4.1}$$

where ¯*r* = *r*/*R*0. Then, according to Equations (8.5) and (8.6), the normalized lithium flux at the particle center is

$$
\vec{J}(0,t) = 0.\tag{A.4.2}
$$

Therefore, the two boundary conditions (8.5) and (8.6) ensure the physical requirement that the lithium flux at the particle center vanishes.

#### **Schriftenreihe des Instituts für Angewandte Materialien**

#### **ISSN 2192-9963**


Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.




Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.







#### KARLSRUHER INSTITUT FÜR TECHNOLOGIE (KIT) SCHRIFTENREIHE DES INSTITUTS FÜR ANGEWANDTE MATERIALIEN

Sodium-ion batteries become an attractive alternative to lithium-ion batteries. Most storage materials exhibit phase changes during operation. The respective phases of a storage material possess different lattice constants giving rise to a strain mismatch which, in turn, causes mechanical stresses and, thus, leads to damage of the electrode particles.

In this work, for non-linear Cahn-Hilliard type models describing diffusion, a thermodynamical framework for the coupling with mechanics is presented. First, a phase-field model for the cathode material NaxFePO4 is studied for the first time to understand phase changes, mechanical deformation, and stress evolution in NaxFePO4. Furthermore, we study the particle size and average concentration dependent miscibility gap of the nanoscale insertion materials LixMn2O4, LixFePO4, and NaxFePO4. Finally, we introduce the nonlocal species concentration theory from a nonlocal free energy density, compare this theory to the Cahn-Hilliard theory, and show how the nonlocality influences the results.

ISSN 2192-9963 ISBN 978-3-7315-1002-4 Phase-field modeling of chemomechanical behavior in battery particles